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ABSTRACT 


Wall functions are often employed to model turbulent flow near solid walls. A 
method has not been available, however, for the application of wall functions to 
generalized curvilinear coordinate systems, particularly those with nonorthogonal 
grids. A general method for this application is developed herein. 

A k — e turbulence model suitable for compressible flow, including the new wall 
function formulation, has been incorporated into an existing compressible Reynolds- 
averaged Navier-Stokes code, F3D. The low-Reynolds-number k-e model of Chien 
(1982) was added for comparison with the present method. A number of features 
were also added to F3D, including improved far-field boundary conditions and viscous 
terms in the streamwise direction. 

A series of computations of increasing complexity was run to test the effectiveness 
of the new formulation. Flow over a flat plate was computed using both orthogonal 
and nonorthogonal grids, and the friction coefficients and velocity profiles compared 
with a semi-empirical equation. Flow over a body of revolution at zero angle of attack 
was then computed to test the method’s ability to handle flow over a curved surface. 
Friction coefficients and velocity profiles were compared to test data. The same case 
was also computed using the Chien (1982) low-Reynolds-number k - e model and the 
Baldwin-Lomax (1978) algebraic model for comparison. All three models gave good 
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results on a relatively fine grid, but only the wall function formulation was effective 
with coarser grids. Finally, in order to demonstrate the method’s ability to handle 
complex flowfields, separated flow over a prolate spheroid at angle of attack was 
computed, and results were compared to test data. The results were also compared 
to the computation of Kim and Patel (1991), in which a, k — e model with a one- 
equation model patched in at the wall was employed. Both models gave reasonable 
solutions, but they require improvement for accurate prediction of friction coefficients 
in the separated regions. 
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1. INTRODUCTION 


1.1 Problem Description 

The understanding of turbulence is of critical importance for the prediction of 
flows encountered in many important engineering applications such as flow over flight 
vehicles, impingment cooling in industrial processes, and the transport of atmospheric 
pollutants. In principle, these flowfields could be predicted by solving the full Navier- 
Stokes equations. This approach is not practical, however, since present computers 
do not have the speed and memory required to resolve the wide range of length and 
time scales in most turbulent flows. In practice, the Navier-Stokes equations are 
employed to resolve large scales, and turbulence models are relied upon to simulate 
the effects of the small-scale motion. 

Turbulence is diffusive, and most approaches to turbulence modeling are directed 
toward computing the rates of turbulent diffusion of momentum and energy. Unfor- 
tunately, a general method for determining these diffusion rates has proven elusive. 
Turbulence models have been developed which work well for certain classes of flows, 
but their range of applicability is limited. Some models, for example, work well for 
attached flows, but perform poorly in regions of separated flow. 

Aside from the generality of turbulence models, another concern is the amount 
of computing power required to apply them. Computations of complex flows may 
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require millions of grid points and hundreds of hours of CPU time, even on the 
fastest available computers. It is therefore important to consider both accuracy and 
computing requirements in the development and application of turbulence models. 

1.2 Historical Review 


1.2.1 Turbulence modeling 

The earliest attempt to analyze the turbulence problem is usually attributed to 
Reynolds (1895). He was trying to explain the result of his famous transition exper- 
iment in which he showed that pipe flow becomes turbulent at a distinct Reynolds 
number. Being familiar with the kinetic theory of gases, Reynolds tried an analo- 
gous approach for fluid flow, decomposing velocities into mean and fluctuating parts. 
When expressions for the decomposed velocities were substituted into the Navier- 
Stokes equations, a set of additional terms appeared. These terms are the gremlins 
which we now call the Reynolds stresses, and the subsequent ninety years or so have 
been littered with attempts to find a general method of predicting their values. 

Since viscous stress in a Newtonian fluid is a linear function of the velocity 
gradient, it was hypothesized that Reynolds stresses behave in the same manner. 
Unfortunately, determining the proportionality constant, the turbulent viscosity, at 
first proved to be as intractable as determining the Reynolds stresses themselves. 

In the 1920s, it was shown that transport equations could be written for moments 
of arbitrary order (Monin and Yaglom 1987). However, each equation for a specified 
moment contains the next higher moment as an unknown. For example, the equations 
for the Reynolds stresses, which are second order moments (the correlation between 
two velocity components), contain third order moments (the correlation between 
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three velocity components) as unknowns. This is the “closure problem” and was a 
harbinger of difficulties to come. 

Some headway was achieved by Prandtl’s “mixing length” hypothesis. It is 
interesting to note that Prandtl, like Reynolds before him, turned toward the kinetic 
theory of gases for inspiration. According to the kinetic theory, kinematic viscosity 
is proportional to the product of a velocity scale (the rms velocity of the molecules) 
and a length scale (the mean free path of the molecules) (Hinze 1987). Treating 
“lumps of fluid” like molecules, Prandtl hypothesized that the turbulent viscosity is 
also proportional to the product of a velocity scale and a length scale. Unfortunately, 
the analogy with molecular motion is on shaky ground at best. Molecules retain their 
identity, while lumps of fluid do not. Also, the length scale of molecular motion is 
small compared to the overall system, and this is not the case for turbulent fluid 
flow (Tennekes and Lumley 1972). Even with these weaknesses, the mixing length 
theory has proven to be useful for the prediction of simple flowfields such as free jets 
and boundary layers on flat plates. Its main drawback is that the proportionality 
constant must be determined empirically, and a given constant is useful only for a 
very limited class of flows. 

An approach very different from mixing length theory was taken by G. I. Taylor 
(1935). Since the Reynolds stresses are expressed as correlations between fluctuating 
components of velocity, it was natural to apply statistical methods to attempt to 
find general expressions for these correlations. Taylor developed this method for 
isotropic and, to a lesser degree, homogeneous turbulence. A great deal of insight 
into the mechanisms of turbulent energy transfer has been gleaned from this work. 
Its application to useful turbulence models has been limited, though, since turbulence 
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is not actually isotropic, and only approximates homogeneity for certain very simple 
flows, such as wind tunnel turbulence behind a grid. 

While statistical methods were being developed, other approaches to improving 
upon mixing length theory were investigated. One of the disadvantages of mixing 
length models is that they do not account for “history” (transport) effects on the 
turbulence. To alleviate this shortcoming, one or more transport equations can be 
employed. It is possible to derive an exact equation for the transport of turbulent 
kinetic energy, although additional unknowns are introduced in the process. The new 
unknowns can be modeled, and the resulting equation can be used to deduce a ve- 
locity scale distribution of the turbulence. Specification of a length scale distribution 
then closes the problem. If the length scale is calculated algebraically, the resulting 
model is known as a “one-equation model,” since one partial differential equation is 
employed. 

One-equation models yield better results than mixing length models for flows in 
which convection and diffusion of turbulent kinetic energy are important (Launder 
et al. 1972). For many complex flows, however, algebraic specification of the length 
scale can be difficult. The next logical step would therefore be to develop a transport 
equation for length scale, or a quantity which can be easily related to a length scale. 
This equation, along with the turbulent kinetic energy transport equation, yields a 
two-equation model. The second equation is usually written for the rate of dissipation 
of turbulent kinetic energy, e, although other quantities are sometimes used, such as 
the length scale, L, (Rodi and Spalding 1970), the rate of dissipation per unit energy, 
u, (Wilcox 1988), and the time scale, r (Abid, Speziale, and Thangam 1991). Two- 
equation models came to the forefront upon publication of a series of papers from Los 
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Alamos Scientific Laboratory (Harlow and Nakayama 1967; Harlow and Nakayama 
1968; Daly and Harlow 1970). Derivation of the second equation is not as rigorous 
as that of the turbulent kinetic energy equation, and this is often cited as a point of 
weakness of two-equation models. Even so, calculation of the length scale as part of 
the model has proven to be advantageous for many flowfields. 

Daunted by the prospect of solving the complete second-moment equations 
and searching for a method to improve the performance of two-equation models, 
Rodi (1972) investigated the possibility of simplifying the second-moment equations. 
He developed an algebraic expression for the Reynolds stresses as a function of the 
dependent variables in his two-equation model, and the model is therefore referred 
to as an algebraic Reynolds stress model. Since the new equation is algebraic, little 
computational effort is required above that for the two-equation model. Although 
algebraic Reynolds stress models show promise, they have not exhibited the expected 
improvements over two-equation models (Ferziger 1987). 

Other variations of two-equation models have also been investigated. One weak- 
ness of two-equation models is that a single velocity scale and a single length scale 
are assumed to be sufficient to describe the turbulence. This implies that the energy 
spectrum is similar in different regions of the flowfield, which is not generally true. 
In “multiscale” two-equation models, the energy spectrum is divided into two parts 
(Launder 1979). The first is the production range, which is the region of highest 
energy. The second is the transfer range, where the energy is transferred from large 
scales to small scales. Separate k and e transport equations are written for each range. 
Multiscale models have shown improvements over standard two-equation models for 
flowfields such as flow over a backward-facing step (Kim and Chen 1989) and swirling 
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jets (Ko and Rhode 1990). The results are not consistently better, however, and a 
significant increase in computer power is required due to the addition of two transport 
equations. 

Another variation of two-equation models is the “nonlinear” model. In some 
flowfields, anisotropy of the normal turbulent stresses is important. An example of 
this is the secondary flow observed to occur in turbulent flow through straight rect- 
angular channels. Since the Boussinesq approximation does not admit anisotropy of 
the normal turbulent stresses, it is impossible to predict these secondary flows with 
the standard model. In nonlinear models (Speziale 1987; Yoshizawa 1988; Barton, 
Rubinstein, and Kirtley 1991), the Boussinesq approximation is replaced by a nonlin- 
ear function of the mean strain rate. This method is not restricted to two-equation 
models, but can be applied to other models which utilize the Boussinesq approxima- 
tion (e.g., algebraic models). Initial results from these models look promising, but 
more applications need to be investigated before their value can be fully assessed. 

As mentioned above, the closure problem precludes the solution of the trans- 
port equations for correlations between fluctuating velocity components. Also, these 
equations contain terms such as pressure- velocity correlations, which are generally 
unknown. Chou (1945) made various assumptions about the unknown quantities in 
the second and third moment transport equations in order to close them, creating 
what is now referred to as a Reynolds stress transport model. An advantage of this 
type of model is that the Boussinesq approximation is not employed. Although the 
Boussinesq approximation is effective for many types of flows, it is known to be in- 
accurate for some flowfields such as wall jets. Chou’s model laid fairly dormant for 
many years, because means for solving the equations for general cases were not avail- 
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able. As computers came into prominence and improved in capability, greater efforts 
were put into the development of Reynolds stress models. These models require a 
great deal of computational effort, and they do not presently yield results which are 
generally better than two-equation models. As they are further refined, it is expected 
that they will come into greater use in the future. 

The goal of all techniques discussed so far is to compute the Reynolds stresses. 
The Reynolds stresses represent momentum transfer averaged over a wide range of 
scales. If a flowfield is computed using a very fine grid, large-scale structures can 
be resolved, and only the momentum transfer occuring at smaller scales needs to be 
modeled. Since the required model represents a subset of the full range of scales, 
it can be simpler in form than models which represent the full Reynolds stresses. 
This approach is called “large eddy simulation.” The disadvantage of large eddy 
simulation is the great amount of computer power required to run with such a fine 
grid. This method is therefore presently constrained to relatively simple flowfields. 

In theory, a grid could be constructed which is fine enough to resolve the full 
spectrum of scales encountered in turbulent motion, obviating the need for any turbu- 
lence model at all. This approach, “direct numerical simulation,” has been applied to 
very simple geometries at low turbulent Reynolds numbers (e.g., Rai and Moin 1989). 
Since a doubling of the turbulent Reynolds number requires an order-of-magnitude 
increase in computer capability (Yakhot and Orszag 1986), it will not be possible to 
use direct simulation to solve “real world” problems in the near term future. It has 
been estimated that if a terraflop (10 12 floating point operations per second) machine 
were available, several hundred thousand years of CPU time would still be required 
to compute a direct simulation of flow over an entire aircraft (Peterson et al. 1989). 
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This would prove to be a major annoyance to typical computer system managers, 
and is therefore untenable. Even so, present direct simulation results are valuable for 
studying the detailed structure of turbulence. Quantities which are not measurable 
can be extracted from the simulation results, and this is an excellent way to check 
details of turbulence models. 

1.2.2 Near-wall modeling 

As solid walls are approached, the structure of turbulent flow changes due to the 
increasing importance (and eventual dominance) of viscous effects. Many turbulence 
models have been developed with the assumption that the flow is fully turbulent (i.e., 
far from walls), and they require additional attention in order to model wall regions 
correctly. 

An early near- wall model which has proven quite useful, and often appears today 
in many guises, is that of Van Driest (1956). Van Driest was looking for a way to 
modify the Prandtl mixing length to account for damping of turbulent eddies near 
walls. He noted that in Stokes’ solution for flow over an oscillating flat plate, the 
amplitude of motion falls off exponentially with distance from the plate. This function 
may be interpreted as quantifying the region of viscous influence. Van Driest used a 
similar function to damp the mixing length near walls, since turbulent effects decrease 
as viscous effects increase. 

As more complex turbulence models came into use, new approaches to modeling 
near-wall behavior were required. Most of these near-wall models attempt to approx- 
imate the effects of anisotropy, which are neglected elsewhere in the flowfleld. These 
models are sometimes referred to as “low-Reynolds-number models,” since they come 
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into play in regions of low turbulent Reynolds number. Harlow and Nakayama (1967) 
presented a tentative anisotropy correction to the turbulent kinetic energy in their 
two-equation model, but they showed no results. Daly and Harlow (1970) used a 
“wall-effect tensor” to modify the fluctuating pressure/strain rate correlation term in 
their Reynolds stress model. They showed that this term drove the peak turbulent 
kinetic energy closer to the wall as the Reynolds number increased, which is in accord 
with experimental data. 

A different approach, “wall functions,” was applied by Patankar and Spalding 
(1970). They reasoned that equations describing the structure of turbulent boundary 
layers, e.g., the law of the wall, could be coupled with numerical solution schemes, 
thereby eliminating the need to resolve the turbulent boundary layer in the region 
where anisotropy is important. This technique is limited by the accuracy and range 
of applicability of the equations employed. 

An early two-equation near-wall model which has been quite influential is that 
of Jones and Launder (1972). They interpreted the e transport equation as modeling 
only the isotropic part of the dissipation rate. Using asymptotic analysis, a term was 
added to the turbulent kinetic energy transport equation to account for anisotropy 
of the dissipation rate. Damping functions were employed for several terms in the e 
equation, and an ad-hoc term was added to bring the maximum level of turbulent 
kinetic energy into line with experimental data. 

Development of both wall functions and low-Reynolds-number models has con- 
tinued in parallel. Chieng and Launder (1980) refined the computation of the wall 
shear stress, and their approach has been implemented by many investigators. Their 
method was further generalized for compressible, separated flow by Viegas, Rubesin, 
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and Horstman (1985). Chien (1982) took an approach similar to that of Jones and 
Launder (1972) to create a low- Reynolds- number model which has gained wide accep- 
tance. There has been a great deal of activity in recent years in the development of 
improved low-Reynolds-number models, usually based on asymptotic analysis. A use- 
ful comparison of eight of these models is given by Patel, Rx>di, and Sheuerer (1985), 
where it is concluded that even the best performing models need more development 
if they are to be used with confidence. Awa, Smith, and Singhal (1990) directly 
compared results of wall functions and a common low-Reynolds-number model for 
three two-dimensional flowfields, and found that wall functions gave comparable or 
better results in all three cases. 

The best choice between the two techniques has yet to be conclusively deter- 
mined. Wall functions yield good results for many problems, and they require less 
computer power than low-Reynolds-number models. Low-Reynolds-number models 
have the potential to be more general and to give better results for some flowfields, 
but that potential has yet to be demonstrated. Both approaches will most likely 
continue to be used in the future. 

1.3 Scope of the Present Research 

One disadvantage of wall functions is that they are difficult to apply to complex 
geometries. Early applications generally involved two-dimensional flows over flat 
surfaces such as duct flows, backward-facing steps, and compression corners. Com- 
putation of flow over complex three-dimensional geometries is now commonplace, but 
a method for applying wall functions to these geometries has not been available. In 
the present work, a method has been developed for the application of wall functions 
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to three-dimensional generalized curvilinear coordinates with nonorthogonal grids. 
A high-Reynolds-number k — e turbulence model with the new wall function for- 
mulation has been added to F3D, a Reynolds-averaged compressible Navier-Stokes 
solver. F3D utilizes an implicit, partially flux-split, two-factor approximate factor- 
ization algorithm, and the ke model utilizes an implicit, fully flux-split, three-factor 
approximate factorization algorithm. The Chien (1982) low-Reynolds-number k - e 
model has also been added for comparison with the wall function formulation. F3D 
contains the Baldwin- Lomax (1978) algebraic turbulence model, which was also run 
for comparison with the present method. 

The new wall function technique was applied to a series of test cases. First, flow 
over a flat plate was computed using two different grids, one which is orthogonal 
and one which is skewed at the wall, to test the nonorthogonal grid capabilities of 
the present formulation. For these cases, the computed friction coefficients were 
compared with those from a semi-empirical equation. Velocity profiles were also 
compared with experimental data. 

Flow over a body of revolution at zero angle of attack was then computed to 
show the method’s effectiveness for flow over a curved surface. Friction coefficients 
and velocity profiles were compared with test data. The same case was also computed 
using the Chien (1982) low-Reynolds-number k - e model and the Baldwin-Lomax 
(1978) algebraic model for comparison. Each of the cases was run on three grids 
with different wall spacings, demonstrating the advantage of wall functions for coarse 
grids. 

Finally, flow over a prolate spheroid at angle of attack was computed using the 
wall function formulation, and results were compared with test data. This demon- 
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strated the effectiveness of the wall function formulation for a complex flowfield with 
regions of separated flow. 
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2. CONSERVATION OF MASS, MOMENTUM, AND ENERGY 

2.1 Introduction 

For turbulent flows, it is not possible to solve the equations of motion numerically 
due to the immense computer power which would be required to resolve the wide range 
of length scales. In order to make the problem tractable, the equations are averaged 
in time, introducing additional unknowns. The additional unknowns, which represent 
turbulent transport of momentum and energy, are then modeled using a combination 
of analysis and empiricism. In this chapter, the technique for averaging fluctuating 
quantities is presented, and it is then applied to the equations of motion. 

2.2 Instantaneous Equations 

The working fluid is assumed to be a homogeneous continuum, and therefore 
may not contain voids or particulates. It is also assumed that the fluid is Newtonian, 
i.e. that the stress is proportional to the rate of strain. Stokes’ hypothesis that the 
bulk and molecular viscosities (A and /z respectively) are related by the equation 
^ = — (2/3)/z is employed. Finally, buoyancy and other body forces are neglected. 
Given these assumptions, the equations of conservation of mass per unit volume, 
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momentum per unit volume, and total enthalpy per unit volume are given by 


dp 5 , . 

at + Tx - {pUi) = 0 

d 0 

fc(P u <) + fa~( pUiU i + 6i > P ~ T '^ = 0 

d d 

fo{pH ~ P) + + 1i ~ “i T ij) = 0 


( 2 . 1 ) 

(2.2) 

(2.3) 


where 


Tij = P 


dui 

dxj 


M: -!*■&) m 

This is a system of five equations with seven unknowns, and must be closed with 
the aid of an equation of state and an expression for molecular viscosity. The fluid 
is assumed to be a perfect gas, 

p = pRT (2.5) 


where temperature is related to total enthalpy by 


H = h + -(«,«,) 


= c p T+-(u iUi ) 


( 2 . 6 ) 

(2.7) 


The molecular viscosity will be calculated from Sutherland’s Law (White 1974), 

p /T\ 3/2 T 0 + S 


P {T> 
Po ~ VT 0 > 


V To / T + S 

where for air, po = 0.1716 mP, To = 491. 6°i?, and S = 199° R. 


(2.8) 


2.3 Averaging Techniques 


Following Reynolds’ approach toward dealing with turbulent flow, values of ve- 
locity and fluid properties are decomposed into mean and fluctuating parts. There are 
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two common techniques of decomposition, “Reynolds averaging” and “mass-weighted 
averaging.” Reynolds averaging is usually employed for incompressible flows, while 
mass-weighted averaging is more convenient for compressible flows. Here, the word 
“average” will be used to refer to averaging over time, defined by 

- 1 f t + At 

* = A th * dT (2 ' 9) 

where <f> is the quantity being averaged, and t is time. In practice, the time must be 
large with respect to the fluctuation time scale, but small with respect to the time 
scale of global changes in the mean flowfield. 

The mass-weighted average is defined by 



(2.10) 


where p represents density, and the tilde is used to indicate a mass- weighted average. 

In Reynolds decomposition, quantites are decomposed into mean and fluctuating 
parts as follows: 

<t> = 4>- \-4> ( 2 . 11 ) 


where <f> is the instantaneous value of the quantity being averaged, <f> is defined by 
equation (2.9), and <ff is the fluctuating part of <f). The analogous equation for mass- 
weighted averaging is 

4> = 4>+<t>" (2.12) 


Now several useful properties of these averages will be developed. First averages 
of fluctuating components will be examined. From equations (2.9) and (2.11), 
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_ 1 rt+At _ 

* = Ttl {4 ’~ ,t ' )dT 

1 rt+ At 1 rt+At_ 

= — / 4>dr — — <t>d,T 

At Jt A tJt 

= <f>- 4> 

= o 

For mass-weighted averages, 

1 rt+At 

P 4/' = — J pTdr 

1 rt+At 

= Ttl M-* 1 iT 

1 rt+At 1 rt+At _ 

= si ^ ~ St i, p4,iT 

= p<f> - p4> 


From equation (2.10), p<f> = p4 > , so 


(2.13) 


p<£" = 0 


(2.14) 


It is important to note that <f/' ^ 0. 

2.4 Mass-Averaged Transport Equations 

The variables in the continuity, momentum, and energy equations will now be de- 
composed into average and fluctuating components, and the results averaged. Equa- 
tion (2.11) will be used to decompose p, p, and r, and equation (2.12) will be used 
for Uj and H. 
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2.4.1 Continuity 


Decomposing the variables in the continuity equation (2.1) and averaging yields 



Applying equations (2.13) and (2.14), 


(2.15) 


d_ 9 . , „ 

S' >+ teJ<^ ) = 0 


(2.16) 


2.4.2 Momentum 


Decomposing and averaging the momentum equation (2.2) yields 



d 


+ 6 ij dT^ p + p,) “ /rr < T 'i + ^ = 0 


dxj 


(2.17) 


Applying equations (2.13) and (2.14), 


+ gH - ( r ij - P<“")| = 0 


(2.18) 


2.4.3 Energy 


Decomposing and averaging the energy equation (2.3) yields 

|[ P (ff + //»)-(?+?')] 

^ _ 

+ «")(# + //") + (^ + q') - (Uj + u'')(T tj + Tjj )] = 0 (2.19) 

Applying equations (2.13) and (2.14), 
d Ft 

Oj(pH ~p) + —[pUjH + pu"H" + q } - UjTjj - u"r 0 ] = 0 


(2.20) 
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It will be useful to 

express H" in terms of h", u„ w", p, and p. 

From the definition 

of total enthalpy, 

rr * 1 

H = h + -UiUi 
£ 

(2.21) 


— h + h! f + — ( U{ + u") (fij + u*l ) 

(2.22) 


= ^ h n + — + uiu*l 

(2.23) 

The definition of mass-average, equation (2.10), applied to the total enthalpy gives 


h = el 

p 

(2.24) 


II 

(2.25) 


~ h ( \p («,+o («.+o 

p 

(2.26) 


- 1 . . 1 pu'lu'! 

- '>+ 2 “.“i + 2 ? 

(2.27) 


Subtracting equation (2.27) from equation (2.23) yields an expression for the fluctu- 
ating component of total enthalpy, 


U" — h" I 1 t/'t/' . - // moo\ 

H = h + -u, u, - - _ + «,«, (2.28) 

Applying this equation to equation (2.20) yields 

| '-(pH - p) + pUjH + pu”h" + q 3 - Ui(T tj - pu'Ju'j) - u'l{r t3 - ^pu'Ju'j) = 0 

(2.29) 

If the boundary-layer approximation is used, the last term on the left hand side may 
be neglected (Cebeci and Smith 1974; Anderson, Tannehill, and Pletcher 1984). It is 
common practice to neglect this term even in full Navier-Stokes computations, and 
this approximation will be used here. 
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The energy equation may be expressed in terms of total internal energy rather 
than total enthalpy. Let the total internal energy per unit mass be denoted by e. 


P 

Multiplying by p, dividing by p, and averaging gives 

pe _ pH p 

P P P 

Using the definition of mass-average and rearranging, 


(2.30) 


(2.31) 


pi = pH -p (2.32) 

Applying equation (2.32) to equation (2.29) and using the approximation 
u i( T ij — \ u '! u 'j) = 0 discussed above, 
f) ft 

dt^ + dx ~ + ijP + pU " }h " + ^ “ P u '! u ")} = 0 (2.33) 

Finally, letting £ denote the total internal energy per unit volume, £ = pe, 
ft ftr/ 

dt^ + d7 j l(^ + p ) + pU> 3 h " + $ ~ “^0 - P u >")] = 0 (2.34) 


2.5 Closure Problem 


The mass-averaged momentum equation (2.18) looks very much like the instan- 
taneous momentum equation (2.2) with the additional term -pu"u”. This term, the 
Reynolds stress, represents the rate of momentum transfer due to turbulent velocity 
fluctuations. Unfortunately, its value is unknown, and the set of equations is no 
longer closed. One’s first inclination might be to derive transport equations for the 
Reynolds stresses. Unfortunately, these equations for second-order moments (aver- 
ages of products of two fluctuating quantities) contain third-order moments, which 
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axe also unknown. Indeed, transport equations for any moments will contain terms 
with higher-order moments. This is the infamous “closure problem.” 

In addition to the Reynolds stresses, there is an additional unknown quantity in 
the energy equation (2.34), pu'-h" . This term represents the transport of energy by 
turbulent velocity fluctuations. 

Both of the unknowns will be computed from a turbulence model, thereby re- 
establishing a closed set of equations. 
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3. k — e MODEL 


3.1 Introduction 


The first task at hand is to model the Reynolds stress term. A coefficient of 
turbulent diffusion may be defined based on the assumption that the turbulent shear 
stresses are proportional to the mean strain rate, in analogy with molecular diffusion. 
An additional term is required to account for the turbulent normal stresses. The 
resulting equation, the Boussinesq approximation, is given by (Anderson, Tannehill, 
and Pletcher 1984) 


— jj-j; ( dui dHj 2. duk\ 2 C ~ . 

~ = "■ (a*J + ~ l 6 ‘ rpk (31) 


where k, the turbulent kinetic energy per unit mass, is defined by 


k = 




(3.2) 


The last term on the right hand side is not included by some authors. Without it, 
however, the turbulent kinetic energy would be identically equal to zero for incom- 
pressible flow. This may be seen by contracting the indices (t = j). 

The unknown Reynolds stress tensor has now been replaced with a function 
of dependent variables from the Navier-Stokes equations as well as two additional 
unknown scalars, the turbulent viscosity [4 t and the turbulent kinetic energy k. 
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It is now assumed that /i t is proportional to the product of a velocity scale and 
a length scale of the large-scale turbulent motion. It is convenient to work with v t , 
which is defined as v t = n t /p, rather than //<. 

u t oc dl (3.3) 


Both scales may vary with space and time. The velocity scale will be taken to be 
equal to the square root of the turbulent kinetic energy, 

t? = (3.4) 


and the length scale will be defined by 


l = 


k 3 ' 2 


(3.5) 


where i is the rate of dissipation of turbulent kinetic energy. Combining equa- 
tions (3.3), (3.4), and (3.5), and defining to be the proportionality constant, 


k 2 

v t = (3.6) 

Assuming that the constant can be determined empirically, knowledge of k and l 
would result in knowledge of the Reynolds stresses. 

There is an additional unknown in the energy equation (2.34), pu"h", which also 
requires some attention. In analogy with the kinematic viscosity, Hinze (1987) defines 
a coefficient of turbulent convective transport for a passive scalar, 


- 7-77 d <t> 

- ^ = v *w, 


(3.7) 


Note that the sign is different on each side of the equation. For mass-averaged 
quantities, it is convenient to define a turbulent diffusion coefficient analogous to the 
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molecular viscosity. For static enthalpy, this coefficient is given by 

77T77 dh 

-p Uj h» = p h — (3.8) 

A turbulent Prandtl number will be defined to relate p h to p t , 

to = -£r (3.9) 

The turbulent Prandtl number is known to vary with space, but a functional re- 
lationship has not been well established (Anderson, Tannehill, and Pletcher 1984). 
Most algebraic turbulence models give good results if the turbulent Prandtl number is 
assumed to be a constant (Anderson, Tannehill, and Pletcher 1984), and this assump- 
tion is also usually made for more complex turbulence models. The commonly-used 
value of Pr t = 0.9 will be adopted here. 

The static enthalpy will now be expressed in terms of the speed of sound a. For 
a perfect gas, 


h 


c p T 


\a 2 

7 R 

a , 2 

7 ~ 1 


Combining equations (3.8), (3.9), and (3.12) yields the final form, 


(3.10) 

(3-11) 

(3.12) 


— pu"h" 


Pt da 2 

Pr t ( 7 - 1) dii 


(3.13) 


The point has now been reached such that the equation set will be closed if the 
k and £ fields are known. Transport equations for k and £ will now be written to 
effect this closure. 
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3.2 Turbulent Kinetic Energy Transport Equation 


The momentum equation (2.2) may be rewritten, replacing the subscript i with k: 


d d 

^P Uk ) + faT.(P UkU i + 6k i p ~ Tk j) = 0 


(3-14) 


Now, multiplying equation (2.2) by Uk, equation (3.14) by adding the results, 
applying the continuity equation (2.1), and simplifying, yields the moment of mo- 
mentum equation: 

d, N . d , N , dp dr jk dp dr jk n , 01KX 

m (pu<u k ) + ^(pu,u,u t ) + u k — - u„— + tti _ - = 0 (3.15) 

This equation will now be averaged. 

£ [p(fi fc + «2) («,- + <)] + P («< + «") (“* + O («i + «")] 

+ <> g - + if = 1 0 (316 > 

Rearranging, and applying equation (2.14), yields 

9 . x d . . dp . . 3p 

a + a ^ {tm,ViUt) + “*a^ + “‘a^ ~ Ut a^ _ “'ai" 


+— (*>«) + (ptiJtiJSj + + pu"u"u t + pu'tu'ju't) 


(3.17) 


The “moment of mean momentum” equation is also needed. This equation may 
be derived in exactly the same way as equation (3.15), but starting with the mean 
momentum equation (2.18), and multiplying by fi* and u, rather than u*. and u t . 



25 


The result is 


d . . d . „ dp _ dp 

at (pUlUk) + + + u ‘aT t 


+ 6 ‘a| (^) +ii ^7 (K<) = o (sis) 


dxj 


Now, subtracting equation (3.18) from (3.17) yields a transport equation for the 
Reynolds stresses: 


h {p«) + 


dxj 


M 9 P 9 P 


— —uv— It" y 4- It" 3 -L 

H dx, Ui dx k +Uk *~ + "< 


' dT tj 
‘ dxj 


dr k} 

dxj 


- pu 'Kjd ~ pu >'’id - 37 - (*<«) 


(3.19) 


: dxj ' 1 ~ J dxj dxj 

Contracting equation (3.19) ( k = i ) and dividing the resulting equation by 2 yields 


d_ 

dt 



(3.20) 


The turbulent kinetic energy per unit mass was defined by equation (3.2). With the 
aid of equation (2.10), its average is given by 


k = (3.21) 

Finally, substituting equation (3.21) into (3.20) yields the equation for the transport 
of turbulent kinetic energy, 
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d_ 

dt 



__3$ 


—u 


+ U\ 


7 


' dxi ' 




tM 

* dxj 



(3.22) 


3.3 Modeled Turbulent Kinetic Energy Equation 


The turbulent kinetic energy equation (3.22) contains unknown correlations on 
the right hand side, so these quantities must be modeled. The first term on the right 
hand side may be expressed as 


„ dp du'-p du'l 


(3.23) 


The second term on the right hand side of the above equation is equal to zero for 
incompressible flow and is expected to be small for compressible flow, so it is common 
practice to neglect it. We now have 


I (>*) + M) 


dxj 


= -£r («&> + K*) + 


dxj 


dlj 


(3.24) 


The first term on. the right hand side of equation (3.24) will be modeled using an 
analogy to equation (3.8). The diffusion coefficient is defined to be p-\- (p. t /(7k), where 
<Jk is a Prandtl number for the diffusion of k. 

d t— . — d \( . nt\ dk' 


d f-r- — d ( n t \ dk 

- — ( u ; P + P u,q = — 1^ + -)^] 


(3.25) 


The third term on the right hand side of equation (3.22) has already effectively 
been modeled by equation (3.1). The final term to be modeled is therefore the second 
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term on the right hand side, the dissipation term. This term will first be split into 


two parts, 


dr^ _d_ ( dvl 

' dxj ax, ^ t. 


(3.26) 


Applying equation (2.4) yields 




" dXj V"* dxj) ' dx j V''“ i dxj 3 6ii d7j 
d ( „a<\ a ( „du’’\ 2 , 

_ du" duj du" diij 2 du" duk 

^ dxj dxj ^ dxj dx{ 3 dx } dxk 

du’; dg d< gg 2 

dxj ^ dxj dxi 3 dij dxk 


d ( „du’j 
+ dZ V UUi d. 


2, d 

3 ° ij dxj 


d ( — jjduk 

% r u ‘ dx k 




(3.27) 


It is generally assumed that compressibility does not affect the dissipation rate. For 
incompressible flow, noting that viscosity fluctuations are equal to zero, the above 


equation reduces to 


ifi g 


d / „d< 

+P a^(“ : &7 


_du" du" _du" du" 
^ dxj dxj ^ dxj dx{ 


(3.28) 


Assuming homogeneity, the first two terms on the right hand side are equal to zero 
(since spacial derivatives of averages of fluctuating quantities are equal to zero). 
Rearranging the remaining terms, and replacing Ji with the equivalent p77 yields 


\dw; 


dxj di{ ) dxj 


(3.29) 


This is the mean density times the dissipation rate, 


(3.30) 
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Finally, applying equations (3.1), (3.25), and (3.30) to equation (3.22) yields the 
modeled turbulent kinetic energy equation: 


d_ 

dt 

+ 


'(g 


_rl diii __ 

ST" 


(3.31) 


The terms on the left hand side represent the total rate of change and rate of convec- 
tion of turbulent kinetic energy per unit volume. On the right hand side are the rate 
of diffusion, rate of production, and rate of dissipation of turbulent kinetic energy 
per unit volume. 


3.4 Modeled Dissipation Rate Equation 

An exact equation for the dissipation rate of turbulent kinetic energy may be 
derived from the momentum equations. This is most commonly carried out by assum- 
ing incompressible flow (e.g. Harlow and Nakayama 1968), although it has also been 
done for compressible flow (El Tahry 1983). An order-of-magnitude analysis of the 
incompressible equation reveals that two terms are of much greater order of magni- 
tude than the others, even though the difference between these two terms is expected 
to be small (Launder 1984). This situation makes it extremely difficult to solve the 
equation numerically. To make matters worse, both terms consist of unmeasurable 
quantities, so even if they could be computed accurately, it would be impossible to 
compare the results with test data. Rather than try to model the exact equation, 
a more heuristic equation is normally used, which mimics the form of the turbulent 
kinetic energy transport equation. This equation, including the compressible terms, 
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is (Coakley 1983) 


a a . a 


(~9c 


+Ci i 


/ a«, 5£t_, 2 

I — + ^ - r<5, 


du k \ 2 


ax, ' ax, z^dxkj ~3 Si i pk 


e 2 


^L-rn- 
dxj 2P k 


C i and C 2 are empirically determined constants. 


(3.32) 
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4. WALL FUNCTIONS 


4 . 1 Background 

The symbol e represents the dissipation rate of turbulent kinetic energy at the 
smallest scales, where the turbulence is nearly isotropic. Near walls, where the tur- 
bulent Reynolds number is low, the turbulence not isotropic, and the dissipation rate 
must be modified accordingly. 

A number of approaches have been used to generalize the model for wall-bounded 
flows. Jones and Launder (1972) created a “low-Reynolds-number” form of the k - e 
model by adding terms to account for the effect of the wall on the dissipation rate. 
This approach has been followed by others, the model of Chien (1982) being notably 
popular. Disadvantages of the low-Reynolds-number models are the additional stiff- 
ness of the equations (Viegas and Rubesin 1983) and the need for fine grid resolution 
near the wall. Also, results from these models are often disappointing when compared 
to experimental data (Chieng and Launder 1980; Patel, Rodi, and Scheuerer 1985; 
Bernard 1986). 

Another approach is to use the high- Reynolds- number version of the k — e model 
away from the wall, and to patch in a different model near the wall. Various models 
may be used in the wall region, such as mixing length models (Lewis and Pletcher 
1986) and one-equation models (Lewis and Pletcher 1986; Rodi 1991). This alleviates 

^rr^ !nvfi(|y PMNy preceding page blank not filmed 
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the stiffness problem encountered with the low-Reynolds-number models, and for 
some flows yields good results. However, since simpler models are used near the wall, 
the concominant disadvantages of these models, such as the need to specify a length 
scale, are encountered. They also require relatively fine grid resolution. 

A third approach, the one to be further developed here, is the use of wall func- 
tions. Wall functions are based on the idea that the basic structure of turbulent 
boundary layers has been well established. Before discussing this structure, some 
definitions are required. In the equations throughout the remainder of the present 
work, all tildes and overbars are dropped except for those indicating correlations be- 
tween fluctuating quantities. An appropriate velocity scale for flow in the near-wall 
region is the friction velocity, defined by 



where r w is the wall shear stress and p w is the density at the wall. Using this velocity 
scale, a nondimensional velocity and a nondimensional length are defined by 


u + = 


u 

i i. 


and 


+ _ u -v 

y 


(4.2) 


(4.3) 


where u is the velocity component parallel to the wall, y is the distance normal to the 
wall, and u is the kinematic viscosity. A typical turbulent boundary layer velocity 
profile, similar to that shown in Anderson, Tannehill, and Pletcher (1984), is shown 
in Figure 4.1. The equation which describes the velocity profile in the log region is 

1 


u 


+ _ 


-ln(y + ) + B 


K 


(4.4) 
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Figure 4.1: Typical turbulent boundary layer velocity profile 


where k is the von Karman constant and B is an additional constant. Values of k 
and B have been empirically determined to fall in the ranges 0.40-0.41 and 4.9-5.5 
respectively (Cebeci and Smith 1974). In the computations in the present study, the 
values of 0.41 and 5.0 are used for k and B. These values were not chosen through 
tuning of results, but simply because they are a frequently chosen pair, and are, in 
the words of Coles and Hirst (1969), “satisfactory and non-controversial.” In the 
viscous sublayer, 

“ + = y + (4.5) 

The equations describing the velocity profile in the inner region are collectively called 
the “law of the wall.” 

Knowledge of the structure shown in Figure 4.1 may be applied in such a way 
that the entire boundary layer need not be resolved numerically. The grid point 
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adjacent to the wall may be placed well away from the wall, and the shear stress 
inferred from the velocity at that point. 

The use of wall functions has several advantages. The other techniques men- 
tioned above require that the entire boundary layer be resolved. The grid point 
adjacent to the wall must therefore be located in the viscous sublayer, typically at 
a y + of less than five. For wall functions, the first grid point is normally located in 
the lower part of the log region, at a y + of approximately 40 to 100. Given that the 
rate at which the grid may stretch away from the wall is limited by most numerical 
solution schemes, wall functions result in a large saving in the number of grid points 
and the amount of computer memory required. Viegas and Rubesin (1983) found 
that approximately half as many grid points were required when using wall functions 
as compared to low-Reynolds-number models. Since the minimum grid spacing is 
much larger for the wall function case, a larger time step may be used for a given 
Courant number (for steady state computations), resulting in further saving in CPU 
time. The reduced memory and CPU required when using wall functions can be 
extremely important for the computation of complex three-dimensional flowfields. 

One disadvantage of wall functions is that the log equation is not accurate for 
some flowfields, such as those with regions of separated flow. Also, the standard wall 
function formulation requires the assumption that the turbulence is in equilibrium 
at the first grid point away from the wall, which is not always the case. Even with 
these limitations, wall functions have been shown to yield results comparable, and 
often superior, to those obtained with low-Reynolds-number models, including com- 
putations of some complex flowfields (Chieng and Launder 1980; Viegas and Rubesin 
1983; Viegas, Rubesin, and Horstman 1985; Chen and Patel 1987; Awa, Smith, and 
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Singhal 1990). 

Given the advantages and disadvantages of each method described above, wall 
functions will be pursued in greater detail. 

4.2 Detailed Formulation 


4.2.1 Introduction 

The first step in applying wall functions is to compute the friction velocity and 
the wall shear stress. The friction velocity is then used to set the boundary conditions 
for k and e at the grid point adjacent to the wall. Finally, the wall shear stress is 
used in the computation of the diffusion term in the Navier-Stokes equations at the 
grid point adjacent to the wall. Development of a general method for applying the 
wall shear stress to the Navier-Stokes equations in generalized curvilinear coordinates 
with nonorthogonal grids is the primary contribution of the present work. 

4.2.2 Friction velocity 


Substituting equations (4.2) and (4.3) into (4.4) and (4.5) gives 



— = -In 1 

U m K ' 

K v J 

(4.6) 

for the log region and 





U 

u*y 



u* 

V 

(4.7) 

for the viscous sublayer. 





Using u and v from the previous time step, the friction velocity u, can be cal- 
culated from equation (4.6) or (4.7). If the grid point adjacent to the wall falls in 
the log region, equation (4.6) is solved using an iterative scheme such as Newton’s 
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method. If the point falls in the viscous sublayer, u m is calculated directly from 
equation (4.7). The appropriate equation is determined as follows. Referring to 
Figure 4.1, equations (4.7) and (4.6) may be seen to intersect at a single value of y+, 
which will be called y+. Neglecting the buffer region as is often done for engineering 
calculations (Tennekes and Lumley 1972), delimits the viscous sublayer and the 
log region. From equations (4.5) and (4.4), 

yt = ) + B ( 4 - 8 ) 

AC 

This equation is solved for y+ using Newton iteration. It is temporarily assumed that 
the point in question is in the viscous sublayer. Using the velocity and viscosity from 
the previous time step, equation (4.7) is solved for tt„, and is calculated from 
equation (4.3). If y + is less than y+, the assumption that the point is in the viscous 
sublayer was correct. Otherwise, the point is actually in the log region, and u, must 
be recomputed using equation (4.6). The wall shear stress may then be computed 
from equation (4.1). 

4.2.3 Boundary conditions for k and e 

The values of k and e require some special attention near the wall. Both quanti- 
ties vary rapidly near the wall, and this region of the flowfield is not resolved due to 
the relatively coarse grid spacing. If k and e can be estimated at the grid point ad- 
jacent to the wall, these values may serve as boundary conditions for the turbulence 
transport equations, and these equations do not have to be integrated all the way to 
the wall. 

For points in the log region, production and dissipation of turbulent kinetic 
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energy are approximately equal, 

P = pe (4.9) 

For a simple two-dimensional boundary layer, letting the y direction be normal to 
the wall, the production term may be written 

du 

P = -pu"v"— (4.10) 

dy 



rearranged, 

_ T t du 
P dy 

From equation (4.6), 

du u. 
dy Ky 


(4.15) 

(4.16) 



Substituting equation (4.16) and the definition of friction velocity (4.1) into equa- 
tion (4.15) yields the final result 


ui 


e = 


Ky 


(4.17) 


If the grid point adjacent to the wall was found to be in the viscous sublayer, 
the assumption that production equals dissipation is not valid. A different method 
must therefore be employed to calculate k and e boundary conditions. The near-wall 
behavior of k may be examined by expanding fluctuating velocity components in 
Taylor series normal to the wall (e.g., Launder 1984). Taking the y direction to be 
normal to the wall, 



(4.18) 


(4.19) 

w " = < + y( 9 y ) w + °(f) 

(4.20) 


where the w subscript refers to the value at the wall. From the no-slip condition, 

= vl = «£ = 0 (4.21) 


Very close to the wall, the flow may be considered incompressible for moderate 
freestream velocities. From the incompressible continuity equation, 



(4.22) 


since (du" Idx)^ and (dw"/dz) w are equal to zero from the no-slip condition. Sub- 
stituting equations (4.21) and (4.22) into (4.18)-(4.20), 
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near the wall. The value of k is known at y = from equation (4.14), and is equal 
to zero at the wall due to the no-slip condition. Equation (4.27) therefore becomes 



(4.28) 


This equation must be used with caution. It has been assumed that the asymp- 
totic analysis is applicable throughout the viscous sublayer, and this assumption is 
questionable. It should be considered an improvement over the assumption that pro- 
duction equals dissipation for points in the viscous sublayer, but is not a definitive 
expression for the turbulent kinetic energy distribution in this region. 

The idea of fitting a parabola to compute k in the viscous sublayer, such as 
equation (4.28), was proposed by Gorski (1986). Rather than fit the parabola up to 
Vc as done above, Gorski fixed the value of the turbulent Reynolds number ( Ret = 
\fkyjv = 20) at the edge of the viscous sublayer, and computed a viscous sublayer 
thickness based on this value. The assumption that Re t = 20 at the edge of the 
viscous sublayer is common practice in finite volume wall function formulations, such 
as that of Chieng and Launder (1980). 
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An expression for e is still needed for points in the viscous sublayer. Following 
Rodi (1991), the length scale equation of Norris and Reynolds (1975) will be used, 

l _ SUL (4.29) 

‘ 1 + 5.3/Re, v ' 


where 


C, = KCj , 3/1 


(4.30) 


Now that k and l t are known, e is simply 


kW 


e = 


(4.31) 


In summary, if the point adjacent to the wall is in the log region, k and e are 
computed from equations (4.14) and (4.17). If it is in the viscous sublayer, k and e 
are computed from equations (4.28) and (4.31). 


4.2.4 Application of r w to the Navier-Stokes equations 

Having computed the friction velocity from equation (4.6) or (4.7), the wall shear 
stress can be computed from equation (4.1) using the density from the previous time 
step. The remaining task is to substitute the wall shear stress into the momentum 
and energy equations. 

The law of the wall was originally developed for two-dimensional boundary layers. 
For three-dimensional boundary layers, if the boundary layer is not highly skewed, 
the wall shear stress calculated above may be divided vectorially into components in 
the two coordinate directions parallel to the wall, proportional to the velocity compo- 
nents. If the boundary layer is skewed, a method of approximating the components 
at the wall, such as extrapolating the velocities, could be employed. Skewness will 
not be considered in the following development. 
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(a) (b) (c) 

Figure 4.2: Definition of 7 coordinate direction 


For generalized curvilinear coordinates, the application of the shear stress de- 
duced from the law of the wall to the momentum and energy equations becomes 
rather complicated. This is due to the fact that the shear stress components in the 
cartesian coordinate system do not necessarily act parallel to the walls. Another 
difficulty is encountered when there is no grid line perpendicular to the wall (i.e. the 
grid is skewed at the wall). The shear stress calculated from the law of the wall 
acts parallel to the wall, in a plane perpendicular to the wall, but this plane is not 
necessarily defined by coordinate directions. A new coordinate direction, called “ 7 ,” 
is defined to be perpendicular to the wall. This is shown in Figure 4.2 for the three 
possible coordinate orientations. 

First, metrics must be established in the new coordinate system. The covariant 
base vectors, which are the base vectors tangent to the coordinate directions, are 
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given by 


a i = 


dr 

W 


(4.32) 


where r is a position vector. These are shown in Figure 4.3. aj represents the base 
vector with components ||, and with similar definitions for the other two 
coordinate directions. Unit covariant base vectors will also be needed, and are given 

by 

(no£) (4.33) 


e, = 


a, 


a, 


The contravariant base vectors are defined by (e.g., Sokolnikoff 1964) 


i _ a 2 x a 3 
ai • a 2 x a 3 
_ 2 _ a 3 x a! 

E — 


(4.34) 


a] • a 2 x a 3 


(4.35) 
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Figure 4.4: Example contravariant base vector 


and 


Q X 3,0 

a = (4- 36 > 

An example contravariant base vector is shown in Figure 4.4. It is sometimes useful 
to express equations (4.34) - (4.36) in vector notation. In foij.C) coordinates, they 
may be represented as 

_ r„ x r c 


V£ = 


Vr? = 


r { ‘ r *? x 

- r C 


r f * r r) x r c 


(4.37) 

(4.38) 


and 


ve = r * x In 
r ( • r„ x r c 


(4.39) 


where r is a position vector, and subscripts indicate partial derivatives. In a compu- 
tational grid, the denominator of equations (4.34) -(4.39) represents the volume of 
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the grid cell, or the reciprocal of the Jacobian of the coordinate transformation, J~ l . 

The 7 direction has been defined to be perpendicular to the wall. As an example, 
refer to Figure 4.2a. In this coordinate system, the 7 direction is perpendicular to 
both the £ and 77 directions. From equation (4.39), it may be seen that VC is also 
perpendicular to the £ and 77 directions. Since the covariant 7 base vector, 


proportionality factor as b, 


/ dx dy 0z\ 
\dj ' d’y' dj) 

(4.40) 

are proportional to one 

another. Denoting the 

dx 

#7 dx ’ 

(4.41) 

dy dC 

d') dy ’ 

(4.42) 

dz h d( * 
07 dz' 

(4.43) 


and 


The 7 direction is now defined, but magnitudes of the 7 metrics still have to be 
established. In the example coordinate system of Figure 4.2a, an equation may be 
written analogous to equation (4.39), 


V7 = -JL01*-. (4.44) 

Tr T n x T i 

Recall that the denominator of this equation is equal to J~ l . Since the values of 
the Jacobians will normally be available in (f , 77, C) coordinates, the Jacobians in the 
(£,77,7) coordinate system will be chosen to be equal to those in the old system. 
This is simply a convenience, so the computation of new Jacobians is not required. 
Comparing equations (4.39) and (4.44), it may be seen that by setting the old and 
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new Jacobians equal, the contravariant base vectors in the 7 and ( directions are also 


equal, 


7x — Cx 
7 y = Cy 


(4.45) 

(4.46) 


7* = Cz 


(4.47) 


J 1 is defined by 


J 1 = ai • a2 x a 3 


or in the coordinate system of the present example, 


(4.48) 


J 1 = • r, x r, 


(4.49) 


Carrying out the vector operations, 


J 3‘t{yr)Z'y + y^(x 7 z n x 7 z 7 ) + z^x 7 y 7 x 7 y^ 


(4.50) 


Substituting equations (4.41), (4.42), and (4.43), and solving for b gives 

, J~ l . 

0 = — — (4.51) 

•^iiy^Cz C y z ri) "i" j/f(Cx 2 r) £t /C z) 4" ^((Xfj^y CxJ^tj) 

The metrics x 7 , y 7 , and z 7 may now be calculated from equations (4.41), (4.42), 

and (4.43). The contravariant base vectors may be calculated from equations (4.34), 

(4.35), and (4.36). In our present example coordinate system, these are 


_ r q x r 7 

r r x r 7 


(4.52) 


Vt,= 


r 7 Xr ( 

r r r„ x r. 


(4.53) 
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The calculation of all the metric quantities that are needed for the coordinate trans- 
formation is now complete. 

Keep in mind that the whole idea of the present procedure is to replace the 
shear stresses in the Navier-Stokes equations with those from the wall functions. For 
finite difference schemes, this means that the shear stresses are required at a point 
between the wall and the point adjacent to the wall (i.e. point 1^). In the absence of 
a streamwise pressure gradient, the momentum equation evaluated at the wall shows 
that the normal gradient of the shear stress is zero at the wall. This means that 
the shear stress at point l| is approximately the same as the shear stress at the 
wall. For cases with streamwise pressure gradients, this is not the case. However, 
it will be assumed here that the wall function shear stress is applicable at point l|. 
This is consistent with the use of the log equation (4.6), which technically is not 
valid for flows with streamwise pressure gradients. It is common practice to use the 
log equation for flows with moderate streamwise pressure gradients, since it gives 
reasonable results for these cases (Launder 1984). Using the pressure gradient from 
the previous time step, a better approximation for the shear stress at point could 
be obtained from the momentum equation if desired. 

The next step is to calculate the shear stresses in physical (x, y, z ) coordinates 
(e.g., r xy ) at point 1^ using the equation 


T 


«7 


(0 + ^) 


Du, duj 
dxj + dii 


2 du k \ 

3 ij dx k ) 


(4.54) 


These should be calculated in the same way that they are calculated in the discretized 
Navier-Stokes equations, i.e. the same averaging procedure should be used to obtain 
values at point l|. The stresses are then transformed to the generalized coordinate 
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system using the standard tensor transformation 


f Q 0 = 


frf fry* r ,- f - 

dx' dxi T 


(4.55) 


(e.g., Sokolnikoff 1964). Here, the "symbol represents quantities in the generalized 
coordinate system. For clarity, Greek letters are used for tensor indices in general- 
ized coordinates, and Roman letters are used for physical coordinates. The tensor 
x represents the physical coordinate x, y, or z, and 7 “ represents the generalized 
coordinate r), or 7 . 

The stresses f a0 do not represent physical quantities. In order to obtain the 
physical components of the shear stresses in the generalized coordinate directions, 
the following equation (Aris 1962) is employed, 


f(ap) = 


QotOt 


(no S ) 


(4.56) 


where f{a(3) represents the physical components of the tensor. g Q0 is the metric 
tensor, 

dx' dx * 


9a0 — 


(4.57) 


$ 7 ° d 7 ^ 

The physical velocity components in the generalized coordinate directions will also 
be needed: 


«(*) = (no 53) 

a 

where U a represents the contravariant velocity components, 


U Q 



(4.58) 


(4.59) 


Capital U is used for the contravariant velocity components and small u is used for the 
physical velocity components in the physical coordinate directions to be consistent 
with standard CFD notation. 
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Figure 4.5: Physical velocity components parallel to wall 

As mentioned above, 3D boundary layer skewness will not be addressed here. 
The shear stresses will therefore be scaled with the physical velocity components. In 
the example coordinate system, the physical velocity component parallel to the wall 
is given by 

V p = V[«(Oe«,+%K] 2 (4-60) 

Here, e €i represents the x, y, and 2 components (corresponding to i = 1,2,3) of 
the covariant unit base vector in the £ direction. Figure 4.5 shows V p for a two- 
dimensional coordinate system (A two-dimensional coordinate system was chosen for 
clarity). Scaling the shear stress components with the velocity components gives 

^7) = -rr^ r M ( 4 - 61 ) 

V P 

nvn) = — ^i 


and 


(4.62) 
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where f| w ] is the wall shear stress from the wall functions. Brackets are used in 
the notation of so that w is not confused with a tensor index. The symbol * 
is used to indicate physical shear stress components which were deduced from the 
wall function equations, as opposed to those computed from equation (4.56). The 
physical shear stress tensor computed from equation (4.56) is now modified to reflect 
the wall function values. In the present example, ?(£y) and f(ryy) calculated in 
equation (4.56) are replaced by the values from equations (4.61) and (4.62). 

The above procedure will now be reversed to transform the new physical stresses 
in the generalized coordinate system to physical stresses in physical coordinates. 
First, an equation is required to calculate tensor components from physical compo- 
nents, i.e. the inverse of equation (4.56). It is convenient here to work in matrix 
notation rather than with tensors. In order to express equation (4.56) in matrix 
notation, matrices will be defined as follows: 



Till) 

r(12) 


p = 

r(21) 

r( 22) 

T (‘ 


_ r(31) 

r(32) 




T 11 

t 12 

T 13 

T 


T 21 

T 22 

r 23 



T 31 

T 32 

r 33 



011 

012 

013 

G 

= 

021 

022 

023 



031 

032 

033 


(4.63) 


(4.64) 


(4.65) 
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y/9n 0 0 

o 0 (4.66) 

0 0 y/g^ 

The ” notation has been dropped, since these definitions hold for any coordinate 
system. Noting that the matrix G is symmetric, equation (4.56) may be written in 
matrix notation as 

P = ST(s- 1 g) T (4.67) 

Solving for T, 

T = S- 1 p[(S- 1 G) T ]" 1 (4.68) 

In the present example the new stress tensor may be computed from equation (4.68), 

t = S- 1 ^>[(s- 1 G) T ]" 1 (4.69) 



where P is 

f(ll) f(12) f(13) 

P= f(21) f(22) f(23) (4.70) 

f(31) f(32) f(33) 

Note that the stress components from equations (4.61) and (4.62) have been substi- 
tuted into the P matrix. Finally, returning to tensor notation (T = f Q ^), the new 
stress tensor must be transformed to physical coordinates, 


-ij dx' dx 1 - Q 8 
t = t 


t j is the new shear stress tensor which is applied to the Navier-Stokes equations. A 


summary of the procedure described above is shown in Table 4.1. 


The application of the new stress tensor to the Navier-Stokes equations is straight- 
forward. The central difference operator 6<t> = <t>j+i/ 2 ~ 4>j-i/2 is typically used for the 
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Table 4.1: Summary of shear stress transformations 


Step 

Actlon Variable 

Equation 

1 

2 

3 

Compute shear stresses in physical coordinates T *j 

Transform stresses to generalized coordinates fa/3 

Compute physical components of stresses in 
generalized coordinate directions ^( a P) 

(4.54) 

(4.55) 

(4.55) 

4 

Split wall shear stress into components parallel 

to the wall *(£ 7 ), f(rry) 

(4.61), (4.62) 

5 

Substitute new shear stress components into : 

physical stress matrix P 

(4.70) 

6 

Compute stresses in generalized coordinates : 

from physical components T 

(4.69) 

7 

Compute new stresses in physical coordinates ^ 0 " 

(4.71) 


diffusion terms. Letting the j direction be normal to the wall, consider the diffusion 
term momentum equation (2.18). Here, r xy represents the total viscous 

and turbulent shear stress. For unit grid spacing, the discretized diffusion term is 

dr xy 

Qy ~ T *yj+ 1/2 ~ T *yj-iii (4.72) 

At the point adjacent to the wall (i.e. j = 2), the value of r xy from the wall function 
calculation is substituted for r Iy ._ 1/a in the above equation. Other stress terms are 
handled in a similar manner. When the equations are transformed to generalized 
coordinates, the shear stresses are still expressed in physical coordinates, so the same 

method applies. For example, the term analogous to equation (4.72) in generalized 
coordinates is 

d 

~ ( J "V*v ). +1/2 - {j-'tyTry) ,_ i/2 (4.73) 

In this equation, r xy from the wall function calculation is substituted for r x iust 

x Vj — 1/2 ’ j 

as in equation (4.72). 
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5. OTHER TURBULENCE MODELS 


5.1 Introduction 

In computing wall-bounded turbulent flows, there are many alternatives to the 
k — e model with wall functions. Two common approaches are low-Reynolds-number 
k — e models and algebraic models. It is therefore desirable to compare their perfor- 
mance to the present wall function formulation. The low-Reynolds-number model of 
Chien (1982) and the algebraic model of Baldwin and Lomax (1978) are quite popu- 
lar and well-tested, and have therefore been chosen for comparison with the present 
formulation. 


5.2 Chien Low-Reynolds-Number Model 

In the standard high-Reynolds-number k — e model, dissipation of turbulent 
kinetic energy is assumed to occur at small scales where the turbulence is nearly 
isotropic. This is reasonable for free shear flows, but is not the case in the vicinity 
of walls. In low-Reynolds-number models, terms are typically added to the standard 
k and e transport equations (Jones and Launder 1972; Chien 1982; Patel, Rodi, and 
Scheuerer 1984; Mansour, Kim, and Moin 1989; Shih and Mansour 1990; Shih and 
Hsu 1991; Michelassi and Shih 1991). The magnitudes of these additional terms drop 
off rapidly away from the wall, resulting in the standard model. 


WTTKTfnNim 


PRECEDING PAGE BLANK NOT FILMED 
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In the Chien model, the dissipation rate is divided into an isotropic part and 
an anisotropic part. The e which appears in the dissipation rate transport equation 
is considered to be the isotropic part of the dissipation rate only. Using asymptotic 
analysis, an expression for dissipation rate near the wall is determined, which is 
intended to account for anisotropy. This expression is then added to the dissipation 


rate appearing in the k transport equation. Also, damping functions are applied to 
the turbulent viscosity and the “destruction of dissipation” term in the e transport 


equation. The resulting equations are 


and 


where 


| (/*) + J- («**) = [(/■ + £) J|] 


+Re 


-l 




2 u duk 

~r k dT t - <* 


{ du{ duj 2 du k \ dui 

\dxj + dx { 3 ,j dxk) Ox, 


—Re 


— 1 2 vpk 

n2 


d , , a 
at i/K) + 


+Re~ 


1 r i\ n (faj , du i 2 C gutM dui 

l k 1 ydxj dx, 3 dx k J j dxj 


2 du k e 2 

-3 C ' 6i ‘ f ’'dT t - C2fp k 


-Re~ l ^e- C * n+ 


n 


u*n 

v 


/= i 




vt 


= (l | -e- c »* I ) 


(5.1) 


(5.2) 


(5.3) 

(5.4) 

(5.5) 


Here, n is the normal distance from the wall. The terms in boxes are those which do 


not appear in the standard high-Reynolds-number k — e model. 
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5.3 Baldwin-Lomax Algebraic Model 

In the Baldwin-Lomax model, the velocity profiles are divided into an inner 
region and an outer region. Each region has its own algebraic expression for the 
turbulent viscosity. The model is given as follows. 


Mt = \ 


(5.6) 


(Mt)inner " — "crossover 
(Mt)outer " • > "crossover 
where n is the normal distance from the wall, and ncrossover is the point at which 
the inner and outer turbulent viscosities are equal. 

The inner formulation is given by 


(Mi) 


inner 


= pl 2 \u 


(5.7) 


where u is the vorticity. The length scale Z is 


l = Kn\l- e (- n+ /' 4+ )l 


(5.8) 


where A+ = 26 and n + is defined by equation (5.3). 
The outer formulation is given by 


) outer — ^ c ^^^wake^kleb 


(5.9) 


where the Clauser constant K c = 0.02688, and = 1.6. Before defining F wai ^ e , an 
additional function is required, 


F = 


zm 

K 


(5.10) 


Tmax is defined as the maximum value of F in the profile, and nmax is the normal 
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distance from the wall at which F m ax occurs, F wa k e is given by 


F wa ke = minimum of < 


(5.11) 


n maxFmax 

Cwfc^maxC^dif) 2 / ^niax 
«dif is the difference between the minimum and the maximum magnitude of velocity 
in the profile, and C w k = 0.25. 

Finally, is the Klebanoff intermittancy factor, 


F kleb = 


1 + 5.5 


( C kleb n N \ 

V n max / 


i -i 


(5.12) 


where = 0.3. 
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6. NUMERICAL METHOD 


6.1 Nondimensional Equations 


For purposes of nondimensionalization, the following reference quantities are 
used: velocity, 0^; density, p temperature, 7^; length, L re /; time, L^f/a^ 
and viscosity, /x^. Using boldface type for dimensional quantities, variables in the 
Navier-Stokes and k — e equations are nondimensionalized as follows: 


= 

OO 


P =P 


h = -jK 


t = 


•Lre// a c 


u 

n- P 


u - Ooo 

P Pocttoo" 

e - 

«oo 

. _ £ 

H-H 

„ X 

Poottoo 3 

aZ 7 


T- T 

T ~iTZ 

T - r 

a- 9 

HooQ'oc / I^ref 

Moo^oo f^ref 


k = 


e = 


a^jLZf 


Applying these definitions to the continuity, momentum, and energy equations 
(2.16), (2.18), and (2.34), and retaining average symbols only in terms containing 
fluctuating quantities, 

|p+£- ( ^) = ° (61) 

(pUi) + [pUiUj + f>ijp - Re~ l (jij - Repu'lu'l)] = 0 (6.2) 

(£ ) + {(£ + p) Uj + Re -1 [qj + Re pu'jh" - u, (r 0 - Repu'/u ")] } = 0 (6.3) 
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where the Reynolds number Re is defined by 

Re =P ooOocLrel 
Moo 

Similarly, the turbulence transport equations (3.31) and (3.32) become 


and 


a Xj 

du{ du 


-\-Re 


-l 


dx 


('*911 

j 2 du k \] dui 2 

j 3 ,j dxk)\ dxj 


duk 

dx k 


P* 


(6.4) 


(6.5) 


d . . d 


J Kg + fe - !«■. £:)]£;- 

The nondimensional turbulent viscosity is given by 

jfc 2 

u t = Rec u — 


dxk 2 k 


( 6 . 6 ) 


(6.7) 


6.2 Vector Form of Equations 


The transport equations for conservation of mass, momentum, and energy are 
usually written in vector form. Let x, t/, and 2 be the three spatial coordinates and 
u, v, and w the respective velocity components. The flux vectors are divided into 
inviscid parts E, F , and G and viscous parts E v , F v , and G v . Equations (6.1), (6.2), 
and (6.3) may then be written 


^ + ^ + ^ + 5*1 D -1 ( dE « , dF « , 9G v \ 

dt + dx* dy* dz \dx + dy + dz ) 


= 0 


(6.8) 
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where the dependent variable vector is 


and the flux vectors are 


E = 


Q = 


P 

pu 

pv 

pw 

£ 


E v = 


“ " 


■ 


-i 

pu 


pv 


pw 

pu 2 + p 


puv 


puw 

puv 

F = 

pv 2 +p 

G = 

pvw 

puw 


pvw 


pw 2 +p 

(£ + p)u 


(£ + p)v 


(£ +p)w 


■ “ 


- 



0 


0 


0 

^cx 


T X y 


T xz 

T X y 

ii 

T yy 

II 

T v* 

^xz 


T v* 


Tzz 

P V 

e 5 


, 




(6.9) 


(6.10) 


( 6 . 11 ) 


where, utilizing the Boussinesq approximation (equation (3.1)) and the model for tur- 
bulent energy diffusion (equations (3.8) and (3.9)), and employing subscript notation 
to indicate partial differentiation, 


[4 2 12 

= (P + Pt) 2 Ux ~ 3 ( v v + w,)j ~ Re-pk 

(6.12) 

= (P + Pt)(Uy + V t ) 

(6.13) 

= (P + Pt) (U z + W x ) 

(6.14) 
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Tyy = in + fit) \v y - ^ (u x + to,)] - Re Ipk 

T yz = (fl + Pt) ( v z + Wy) 

r zz = {p + fit) jw 2 - | («, + »„)] - Re ^pk 

e l = ur - + VT *y + WT ** + [(^TI) (j; + p^r) (° 2 )J 

ft = UT X y + VTyy + WT zy + ^ ^ (a') ^ 

gl = UT XZ + vr yz + WT ZZ + ^T[] ( P r + ~pr) (° 2 ) J 

The turbulence transport equations may be put in a similar form, 


dQ t dE t dF\ dG t _ _ x ( dE* OF* dGtv \ 
dt + dx + dy + dz 6 \ dx dy dz ) 


= H t 



( 6 . 15 ) 

( 6 . 16 ) 

( 6 . 17 ) 

( 6 . 18 ) 

( 6 . 19 ) 

( 6 . 20 ) 

( 6 . 21 ) 

( 6 . 22 ) 

( 6 . 23 ) 

( 6 . 24 ) 

( 6 . 25 ) 

( 6 . 26 ) 
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Before writing the source term H t in vector form, it is useful to define a symbol 
for a term which is proportional to the production of turbulent kinetic energy, V, 

t, 2 , du k 


V = Re' 1 
Now H t may be written as 


( du { 

dui 

2 du k V 

r U; 

+ 

CD 

£ V 
1 



dx k 


H t = 


V — pe 


c^v-c^ 

6.3 Coordinate Transformation 


(6.27) 


(6.28) 


The standard transformation to generalized curvilinear coordinates (e.g. An- 
derson, Tannehill and Pletcher 1984) will be employed here. Transformed time is 
identical with physical time, and transformed spatial coordinates are general func- 
tions of physical spatial coordinates and time, 


T = t 

f = £(x,y,z,t) 

v = 

C = ((x,y,z,t) 


The Jacobian of the transformation is 


djtrpQ 

d(x,y,z) 


(6.29) 

(6.30) 

(6.31) 

(6.32) 


(6.33) 


Partial derivatives of quantities with respect to physical coordinates may be expressed 
in terms of transformed coordinates through the chain rule. Letting (p represent a 
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variable to be differentiated, 


4>t — <f> T + 4- 4- (6.34) 

4*x — d - Vx^r) 4" Cx^c (6.35) 

4>y £j/0{ 4” Vyfir) 4" Cy^C (6.36) 

4>z — + Vz^r) 4- Cz<Ac (6.37) 


Applying the chain rule to equation (6.32), dividing by the Jacobian, and rearranging 
and cancelling terms results in the transformed set of equations. 


dQ dE dF dG _ x (dE v , dF v , 8G V \ 

aT + af + a^ + W~ Re {-W + -^ + w) 


= 0 


(6.38) 


where 


P 

pu 


Q = J~ l 


pv 


pw 

£ 


E = J -1 


P U 

puU + £ x p 
pvU + i y p 
pwU + ( z p 
(£ + p)U - hp 


(6.39) 


(6.40) 
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P V 

puV + TJ x p 

PVV + T]yP 

pw V + r) z p 
{£ + p)U - Tjtp 


pW 

puW + C x p 
pvW + CyP 
pwW + ( z p 

(£ + p)U - C tP 


0 


^xT~xx 4" (,y^xy 4" (,z T xz 
£x T xy 4" Cy T yy 4“ ^z^yz 
£x 7~xz 4" Cy ^~yz 4" ^z 7"zz 


o 

Vx^xx “ 1 “ Vy^zy Vz^xz 
r )x'T X y + TjyTyy + T\zTyz 
T) x^~xz “1“ T}yTyz ”f“ T)z7~zz 
VxCl + T) y f* + 7 i z gl 


(6.41) 


(6.42) 


(6.43) 


(6.44) 
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where 


0 



Cx 7~xx “I" Cy T xy 4" Cz^~xz 
Cz^xy 4” Cy T yy 4" C z T yz 
Cx^xz 4” Cy r yz 4" C z^zz 
( x el + ( y fZ+Cz9Z 


(6.45) 


T xx 


Txz 

T VV 


'yz 

Tzz 


e 


V 

5 


(p 4- fit) 4- VxUr, + Cr«c) 

- \ [(Zv V t 4- T)yV n 4- <>c) 4- (&«>( + VzWr, + <z™<)]} 
2 

-•Re -pk 


(p + Pt) [(£„«{ + T) y U v + CyUc) 4- (ZxVt 4- VxVr, + <><)] 
{fi 4- fit) [(&«{ 4- 9zU v 4- Cz«<) 4- (£*w € 4- Ijxti), + Cx«^c)] 

(M 4 -fit) + TJyVr, + CyV() 

- ^ [(&«( 4- T Ixtlr, + Cx«c) 4- (Cz^ 4- + Cz^c)]} 

2 , 

-Re-pk 


(p + Pt) ((&v* 4- T] z v v + Cz^c) 4- (&«>( + W + Cy^c)] 
(p 4- Ah) ||(6«>* 4- Vz^T, + Cz«>c) 

“ I [(&u ( + »!,«, 4- <,«c) 4- (£„«* 4- 4- »fct; ( )]} 

„ 2 , 

-Re-pk 

ut xx + vr IV -f wt X2 


h rjj (& + £;) [«'(“ 2 )< + -»*(“ 2 ), + 


(6.46) 

(6.47) 

(6.48) 


(6.49) 

(6.50) 


(6.51) 

(6.52) 
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fb ur xy + vr yy ~t" wr zy 

Q§ — UTxz “H V7~yz im'zz 


+ 


(7-1) 

and the contravariant velocities are 


(6.53) 


(6.54) 


U = Z t +£ x U + iyV + t Z W 

(6.55) 

V = lit + r) x u + r)yV + T} z w t 

(6.56) 

W = Ct + Cx u + (y v + Cz w t 

(6.57) 


The transformed turbulence transport equations are 


3Q, 8B, 9F, 9G, 

-W + ST an + IT" 

where the dependent variable vector is 


■( 


dE tv dF tv dGtv 

d£ drj d( 


)■ 


H t 


Qt- J 


-i 


The flux vectors are 


E t = J 


-l 


F t = J 


-l 


G t -J 


-l 


pk 

PC 

C/pfc 

Upe 

Vpk 

Vpe 

Wpk 

Wpe 


(6.58) 


(6.59) 


(6.60) 


(6.61) 


(6.62) 
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Etv — J 


;("+£) [(£+«?+«?)£ 

-1 | + (6>>, + Of, + 6»fc) ^ + (60 + 60 + 60) 

?(" + £) [(G 2 +e?+o 2 )f 

+ ( 61 , + 6 l» + 6 >),) + (60 + 60 + 60) 


F,„ = j 


; (f + S) [(6% + 6>fc + 6%) ^ 

-1 I + (rf + T f V + ^i) ^ + (>liG + TJyCy + 7 7,0) if] 


; (/< + £) [(6>i» + 6% + 6>fc) ^ 

+ (ll + 'f, + r /j) ff + (ifeO + 1,0 + 1,0) j J 

p (f + “) [(60 + 60 + 60 ) 

+ (>?,o + >i,o + i«o> ^ + (a + o 2 + o 2 ) ^ ] 

} (/* + £) [(60 + 60 + 60 ) If 

+ {>iiO + >1,0 + TfcO) + (c, + C 2 + C 2 ) if] J 

and the source term vector is 


<?,„ = J- 1 


H t = J- 1 


V — pe 

C^v-c^ 


(6.63) 


(6.64) 


(6.65) 


(6.66) 


6.4 Navier-Stokes Solver 


Equation (6.38) will first be semi-discretized in time. Using first order forward 
differencing, 


Q n+l -Q n + At 


— Re' 


dE n+1 dF n+1 dG n+1 
dk + drj + ~dc~ 

1 (dE: +1 dF v n+1 dG n v +i \ 

{ dt + drj + d( )\ 


= 0 


(6.67) 



67 


where superscripts refer to time level. The flux terms must now be linearized. Ex- 
panding E n+1 in a Taylor series, 

E n+ 1 = E n + At (^j +0 [(At) 2 ] 

Using the chain rule, 

Substituting equation (6.69) into (6.68) and discretizing (dQ/dij n yields 

E n+ 1 = E n + A n AQ n + O [(At) 2 ] (6.70) 


( 6 . 68 ) 


(6.69) 


where 


A Q n - Q n+1 - Q" 


(6.71) 


and 



(6.72) 


Jacobian matrices may also be defined for the other inviscid fluxes, B n 


( 3F/3Q )" 


and C — (dG/dQ} , and similarly for the viscous fluxes. Analytical expressions for 
the inviscid and viscous flux Jacobians are given in Appendix A. 

Applying equation (6.70) to equation (6.67) results in the delta form of the 
equation, 


{/ + At [(^A" + 3,5" + 3 f C") - Re~ l (c^ + + 3 ( C V ")]} AQ 

= -At [( 3 cE n + d^F" + 3 f G") - Re- 1 (d^ + d.F? + 3 C G?)] (6.73) 

where operator notation is now being used for the partial derivatives. Equation (6.73) 
is solved using an implicit, partially flux-split, two-factor approximate factorization 
algorithm (Ying et al. 1986). In the original algorithm, all cross-derivative terms 
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were neglected, and viscous terms were retained only in the direction normal to the 
wall (thin-layer assumption). Here, viscous terms are retained in all three directions, 
with all cross-derivative terms neglected. In the flux-split direction (£)> the viscous 
terms cannot be linearized if a tri-diagonal solver is to be employed, so they are 
included only on the right hand side. 

The flux vector splitting method of Steger and Warming (1981) is used for con- 
vection terms in the £ direction. The £ direction inviscid Jacobian A is therefore split 


as follows, 




A = A + + A~ 

(6.74) 

where 




A + = f A+f- 1 

(6.75) 

and 




A" = f A-f- 1 

(6.76) 


A + is a diagonal matrix containing the positive eigenvalues of A,, and A contains 
the negative eigenvalues. The columns of T are the right eigenvectors of A,- . 

Substituting difference operators, adding smoothing, and applying approximate 
factorization to equation (6.73), 

4- At [V { (•<4 + ) + 6<C n — Re — ^imp<] } 

x {/ + At [a^(A ) + Sr)B n ~ R e &r)B v — T^jniptj]} A Q 
= A t [V f (£ + )” + A 1{ (£-)" + 6,F" + f ( G" 

—Re'* — (pexpn + F’expc) Q n ] (6-77) 

The difference operators are defined by 


= 4>j ~ <j>j - 1 
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= <f > j+1 - 

= ^(0i+i - 

6^4> = <t>j+\/2 — (j)j- 1/2 (6.78) 


Similar definitions hold for the other two coordinate directions. Smoothing is required 
for the central-differenced directions to ensure stability. The smoothing operators are 


given by 


T> ■ = 7 -1 

impc J 

^expc = J 1 


-(- 2.5e 4 6<; j ^ 

+ € ^c $1 


<J 

j 


j 


(6.79) 

(6.80) 


where 



(6.81) 


e 2 and e 4 are the second and fourth order smoothing coefficients, and is the spectral 


radius of the inviscid flux Jacobian matrix. The factor of 2.5 in the implicit smoothing 


operator is included to ensure stability, given the fourth order explicit smoothing. 
The p function is designed to switch from fourth order to second order smoothing in 
regions of rapid spatial variation in pressure, such as at a shock wave. The spectral 
radius scaling is included to emulate the numerical dissipation inherent in upwind 
difference schemes (Pulliam 1985). The rj direction smoothing operators are similar 
to equations (6.79) and (6.80), with all occurrences of the symbol £ replaced by rj. 


6.5 k — e Solver 

The k-e equations are solved uncoupled from the Navier-Stokes equations. The 
linearized equations are simpler in form than the Navier-Stokes equations because 
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they are coupled only through the source terms if the functional dependence of on 
k and e is neglected in the viscous terms. Flux Jacobians for the convection terms 
are simply diagonal matrices containing the contravariant velocities. The linearized 
equations in delta form axe 


{/[i + Af(d{t/ + a„v + a c wo] 

-A tRe~ l (J-WPAJ) + dr, (j- 1 N0 A d 1 ,j) + {j~ x M0A ^)] 

-Atl)^} A Q t = -At [(^£ t n + dr,F? + 3<G?) 

-Re' 1 + 5„F» + d c Gl) - H?] (6-82) 

where 


A7 = 




(6.83) 


and 


(h=e x +e y +g 

(6.84) 

fa = v z + nl + 

(6.85) 


A> = Cx + Cj + Cz 


( 6 . 86 ) 


An analytical expression for D™, the turbulent source term Jacobian, is given in 
Appendix A. 

Equation (6.82) is solved using a conventional three-factor approximate factor- 
ization scheme with flux vector splitting and first-order differences for the convection 
terms in all three directions. Central differences are used for all diffusion terms. 
Several methods of handling the source terms are available. Here, they have been 
placed in a single factor, since this method is both simple to implement and is also 



71 


computationally efficient (Shih and Chyu 1991). The resulting equation is 

{/ [l + At (V f £/ + + A ( tr)”] - At Re~‘ [i s j)] } 

x {/ l 1 + At (V,V + + A,y)"] - At Re~' |J, {j-'MM'j))} 
x {/ [l + A( (v c W* + A C I V-)"] - At Re-' p< {j-'MPeh ./)] 

- AtT»"}A4 

= -A({[V { (£+)" + A ( {E 7 )" + V,(F/)" + A, (/•,-)" 

+V < (6', + )" + A C (G'-)"] -Re-' (? { £J + b„F" + « c G;') - H?} (6.87) 

The ^ and tj direction factors require solution of uncoupled equations, since there 
is no source term Jacobian present. A specialized solver was therefore written for 
banded tridiagonal matrices to enhance computation speed. In the ( direction, a 
block solver is required due to the source terms. Since two-by-two blocks may be 
inverted algebraically, a second specialized solver was written for the £ direction. 
Details of both of these solvers are given in Appendix B. 

An additional consideration in the solution of these equations is the possibility 
of obtaining negative values for k and e. Negative values are physically impossible, 
but are admitted by the modeled transport equations. Lower limiters are therefore 
used for both equations. After each step, any values of pk or pe that are below the 
specified limit are bumped up to that limit. The freestream values of pk and pe are 
used for the limiters, since the physical (as opposed to numerical) values are not 
expected to fall below those in the freestream. One run was started from scratch (all 
dependent variables set to freestream values) with the limiters turned off, and it was 
unstable. Addition of the limiters solved the problem. After the solution had partly 
converged, the limiters were again turned off, and the solution remained stable. The 
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starting transient had caused the unphysical values in the first case, and once the 
solution settled down, the limiters were unnecessary. 
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7. RESULTS 


7.1 Introduction 

A formulation designed for general geometries and skewed grids must also work 
on simple geometries and orthogonal grids. The ubiquitous flat plate has therefore 
been chosen for the first set of test cases. First, turbulent flow over a semi-infinite 
flate plate was computed on an orthogonal grid. The next step was to verify that the 
formulation is effective when the grid is skewed at the wall, so the same flat plate 
was solved with a skewed grid. 

The tensor transformations in the present formulation contain functions of all 
of the metrics. When any of the generalized coordinate directions are orthogonal to 
physical coordinate directions, some of the metrics will be identically equal to zero. 
One final flat plate test case was run with the entire domain at a thirty degree angle 
with the physical coordinate system (at zero angle of attack) to bring additional 
metrics into play and further test the method. 

One additional flat plate test case was run to verify the coding of the Chien 
(1982) low-Reynolds-number model. 

After having verified the present formulation for flat surfaces, the next step 
was to solve a well-behaved (non-separated) flow over a curved surface. Ramaprian, 
Patel, and Choi (1978, 1980) took measurements of the turbulent flowfield around a 
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body of revolution at both zero and fifteen degree angles of attack. The zero degree 
(axisymmetric) flowfield was computed using the present wall function formulation as 
well as the Chien (1982) low-Reynolds-number model and the Baldwin-Lomax (1978) 
algebraic model, and the results compared to the measurements. Each computation 
was carried out using grids with three different wall spacings to determine the grid 
sensitivity of each model. 

Finally, the present wall function formulation was tested on a complex three- 
dimensional flowfield. Separated flow over a prolate spheriod at ten degrees angle 
of attack was computed, and the results compared to the measurements of Kreplin, 
Vollmers, and Meier (1982). 


7.2 Flat plate 

The first test case is the computation of incompressible, turbulent flow over a 
semi-infinite flat plate with zero pressure gradient. The freestream Mach number 
was set to 0.2 to minimize compressibility effects without getting too close to the 
incompressible limit of the solver. The Reynolds number based on freestream velocity 
was 1 x 10 6 at the upstream boundary, and 8 x 10 6 at the outflow boundary. The 
reference length was defined such that the distance from the virtual origin of the 
boundary layer to the upstream boundary was one unit of length. All lengths were 
normalized by this distance. 

The streamwise and normal coordinate directions are x and z in physical coor- 
dinates and 77 and ( in transformed coordinates. Wall functions are expected to give 
the most accurate results when the grid point adjacent to the wall falls in the log 
region. The distance from the wall to the grid point adjacent to the wall, which will 
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be referred to as the “wall spacing,” was chosen to equal 0.0035 (nondimensionalized 
by the reference length described above) so that A z+ « 140, which is well into the 
log region. A is the wall spacing in wall coordinates. The domain extended from 1 
to 8 in the streamwise direction, and 0 to 1.5 in the normal direction. The location of 
the outer edge was chosen to exceed ten times the estimated boundary layer thickness 
at the downstream boundary. The 101 x 61 (streamwise x normal) orthogonal grid 
is shown in Figure 7.1. 

An upstream streamwise velocity profile was specified using the law of the wake 
(Coles 1956), 

u + = -ln{y + ) + B + -W (7.1) 

K K> \0 / 

where 

M '(!H sin2 (i!) (7 - 2) 

II is a function of the streamwise pressure gradient, and is approximately equal to 
0.5 for pressure gradients of zero. The boundary layer thickness, <5, was estimated 
from the equation (White 1974) 

- « 0.37fleJ 1/5 (7.3) 

X 

At the inflow boundary, the normal velocity component was set to zero, the den- 
sity was fixed at the freestream value, and the pressures were set by extrapolation 
(p^j = p v - 2 ). At the outflow boundary, pressure was fixed at the freestream value, 
and density and both components of momentum were extrapolated. At the outer 
edge, density, both components of momentum, and pressure were extrapolated. At 
the wall, velocities were set to zero, and density and pressure were extrapolated. 

The upstream turbulence quantities were more difficult to estimate. A detailed 



76 



(a) Complete grid 
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turbulent kinetic energy profile was measured for incompressible flow over a flat 
plate with zero streamwise pressure gradient by Klebanoff (1955), but at a Reynolds 

number of 4.2 x 10 6 . These data were used to estimate the distribution at the present 
Re = 1 x 10 6 . 

The friction coefficient was first calculated using the equation (White 1974) 


C f « 


0.455 

In 2 (0.06Re z ) 


(7.4) 


Combining the definition of friction velocity, equation (4.1), and the definition of 
friction coefficient, 


1 

the friction velocity can be expressed as a function of friction coefficient, 


(7.5) 


u , 




(7.6) 


Since the grid spacing at the wall was chosen such that the first grid point away from 
the wall was in the log region, k at that point was calculated from equation (4.14). 

The next step was to interpolate the Klebanoff k distribution onto the present 
grid using a cubic spline (Hornbeck 1975). As was expected, the value of A: at the point 
adjacent to the wall did not match the value calculated from equation (4.14). The 
whole k distribution was multiplied by the ratio of the two values to force the correct 
value at the given point. Finally, freestream values were cut off at k/U^ = 0.0002, 
Klebanoff’s approximate freestream k. 

The upstream e distribution was estimated using a method described by Launder 
et al. (1972). The ratio of turbulent shear stress to turbulent kinetic energy was 
estimated to equal the constant 0.3. Since the k distribution was estimated above, 
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the turbulent shear stress distribution was also “known.” Finally, the e distribution 
was computed from equation (4.13). 

At the outflow boundary and the outer edge, k and e were extrapolated. Values 
of it or e are not required at the wall when using wall functions. 

The friction velocity was computed at each point along the plate from either 
equation (4.6) or (4.7). Since the friction coefficient is easily deduced from the fric- 
tion velocity, comparison with experimental friction coefficients is a good check on 
the effectiveness of a wall function formulation. The friction coefficient distribution 
for incompressible turbulent flow over a flat plate with zero pressure gradient is well 
established. An accurate and well-tested equation for the friction coefficient distri- 
bution is given by (White 1974) 


Re x = -^A 4 + - 


kB 


12 


K 


e*(Z’-4Z + 6)-6-2Z-^-§^ 


(7.7) 


where 

A = v /2 ]C, 


(7.8) 


and 


Z = k A 


(7.9) 


The computed results are within 1 of the values from the equation. The 
computed friction coefficients are compared with equation (7.7) in Figure 7.2. The 
computed results are within l|% of the values from the equation. 

A sample velocity profile is compared with test data in Figure 7.3, where 6 is the 
boundary layer thickness and Uinf is the freestream velocity. The measurements were 
taken by Klebanoff (1955) at a Reynolds number of 4.2 x 10 6 based on freestream 
velocity and distance from the virtual origin of the plate. The computed profile 
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Figure 7.2: Friction coefficient, flat plate, orthogonal grid 



U/Uinf 

Figure 7.3: Velocity profile, flat plate, orthogonal grid 
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is at the same Reynolds number assuming that the virtual origin is at x = 0. The 
boundary layer thickness for the test data was taken to equal the “nominal thickness 
of three inches quoted in Klebanoff’s paper. For the computation, the boundary layer 
thickness was taken at the point at which the velocity reached 99% of the freestream 
velocity. 

At the grid point adjacent to the wall, the computed velocity is slightly higher 
than the test data. Since this grid point falls in the log region, the computed friction 
coefficient must also be larger than the value corresponding to the test data. It 
would be more accurate to compare velocity profiles at the same momentum thickness 
Reynolds number rather than Reynolds number based on distance from the virtual 
origin of the plate, since it is difficult to locate the virtual origin precisely. For the 
present computation, however, it is also difficult to compute the momentum thickness 
accurately due to the coarse grid spacing at the wall. Since the boundary layer is not 
resolved near the wall, the numerical integration required to compute the momentum 
thickness would also be inaccurate. 

7.2.1 Skewed grid 

In the present wall function formulation, a new coordinate direction is defined 
normal to the wall, and physical shear stresses are computed using this new direction. 
A simple test case with a grid that is nonorthogonal at the wall was desired to check 
this part of the procedure. A case similar to that computed above was chosen, 
but using a grid which is skewed near the wall. The grid lines parallel to the plate 
(constant () are identical to those in the orthogonal grid. The grid is orthogonal to the 
plate at the upstream boundary. Moving downstream, the grid lines gradually become 
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more and more skewed, until they reach an angle of sixty degrees with respect to the 
plate. They then remain at this angle for the rest of the domain. This configuration 
was chosen for two reasons. First, the upstream boundary definition was identical to 
that for the orthogonal grid case, and it was therefore easy to implement. Also, some 
of the metrics (e.g. dC,/ dx) are changing along the transformed coordinate directions, 
so additional terms are brought into play in the computation. The grid is shown in 
Figure 7.4. The freestream conditions and boundary conditions are identical to those 
for the orthogonal grid case. 

The resulting friction coefficient distribution is shown in Figure 7.5. It is iden- 
tical to that for the orthogonal grid, verifying that the shear stresses are being com- 
puted correctly in the new coordinate syatem. 

The velocity profile at a Reynolds number of 4.2 x 10 6 is shown in Figure 7.6. 
It, too, is identical to that in the orthogonal grid test case (Figure 7.3). 

7.2.2 Angled domain 

Since the geometry of the flat plate test cases is simple, many metric terms used 
in the transformation of the shear stresses are identically equal to zero (for example 
dy/dx). The flat plate test case with the orthogonal grid was therefore repeated, but 
the entire domain was angled thirty degrees with respect to the physical coordinates. 
The plate was still at zero degrees angle of attack. In this computation, none of 
the physical coordinates were orthogonal to the transformed coordinates, so it was 
a more general test of the wall function formulation. This case was run at the same 
freestream conditions as the previous cases. The resulting friction coefficients are 
shown in Figure 7.7. The angled domain has not affected the results, verifying that 
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(b) Close-up of region of changing skewness 
Figure 7.4: Skewed flat plate grid 



(c) Close-up of downstream boundary 
Figure 7.4 (Continued) 



1e+06 2e+06 3e+06 4e+06 5e+06 6e+06 7e+06 8e+06 

Re* 

Figure 7.5: Friction coefficient, flat plate, skewed grid 
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U/Uinf 

Figure 7.6: Velocity profile, flat plate, skewed grid 

the shear stress transformations axe working correctly for this case. 

The velocity profile at a Reynolds number of 4.2 x 10 6 is shown in Figure 7.8. 
It is identical to those in the previous two test cases. 

7.2.3 Low-Reynolds-number model test case 

In addition to the flat plate test cases for wall functions, a test case was run to 
check the coding of the Chien (1982) low-Reynolds-number model. For this model, a 
fine grid is required at the wall in order to resolve the buffer region and the viscous 
sublayer. In these regions of the flowfield, which are not resolved when using wall 
functions, k and e vary rapidly, and the method used to estimate the upstream k and 
e distribution in the previous test cases becomes questionable. In order to simplify 
specification of the upstream boundary conditions, a semi-infinite flat plate with the 
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Re x 

Figure 7.7: Friction coefficient, flat plate, angled domain 



U/Uinf 


Figure 7.8: Velocity profile, flat plate, angled domain 
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leading edge immersed in the freestream was used. The 91 x 91 grid is shown in 
Figure 7.9. For this grid, the reference length was the distance from the leading edge 
to the downstream boundary. The grid spacing at the wall was 1 x 10 -5 based on the 
reference length. The domain extended from x = —0.5 to x = 1.0 in the streamwise 
direction and z = 0 to z = 0.5 normal to the plate, with the leading edge at x = 0. 
The Reynolds number (based on freestream velocity) was 9 x 10 6 at the downstream 
boundary, and the freestream Mach number was 0.2. 

At the upstream boundary, the freestream density and momentum were fixed and 
the pressure was extrapolated. As in the previous cases, the freestream turbulent 
kinetic energy was fixed at k/U^ = 0.0002. The freestream dissipation rate was 
computed based on an assumed nondimensional freestream turbulent viscosity of 1. 
At the stagnation streamline, all quantities were extrapolated. Specification of the 
wall, edge, and downstream boundary conditions were the same as in the previous 
test cases. 

It is preferable to plot friction coefficient versus Reg (momentum thickness 
Reynolds number) rather than Re x (Reynolds number based on distance from the 
leading edge) in order to avoid difficulties in locating the virtual origin of the plate. 
This was not possible in the previous test cases due to the coarse grids, but was 
appropriate here. A semi-empirical equation for the friction coefficient distribution 
is given by White (1974), 

Re, = (3.75 - (7.10) 

where A is defined in equation (7.8). The computed friction coefficient distribution 
is compared to the values from equation (7.10) in Figure 7.10. Simpson’s rule was 
used for the numerical integrations to calculate the momentum thicknesses. The 
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(b) Close-up of leading edge region 
Figure 7.9: Flat plate grid, low-Reynolds-number model test case 
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Figure 7.9 (Continued) 



Figure 7.10: Friction coefficient, flat plate, low-Reynolds-number model test case 
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U/Uinf 

Figure 7.11: Velocity profile, flat plate, low-Reynolds-number model test case 


agreement is good except near the the leading edge, where the boundary layer is too 
thin to resolve. 

The boundary layer profile is compared to the Klebanoff (1955) data in Figure 7.11. 
The momentum thickness Reynolds number of the Klebanoff profile was calculated, 
and a profile was chosen from the computation at approximately the same momen- 
tum thickness Reynolds number. Figure 7.11 shows that the computed velocities are 
somewhat low near the “corner” of the profile, and that the match is good toward 
the outer edge. At first glance it appears that the momentum thicknesses cannot be 
equal for both curves. All of the computations exhibited a slight overshoot at the 
edge of the boundary layer, so this is probably the cause of the discrepancy. 
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7.3 Body of Revolution 


7.3.1 Introduction 

The next step was to test the wall function formulation on a curved surface. A 
well-behaved (i.e., non-separated) flowfield was desirable, since the purpose was to 
check the stress tensor transformations, and not yet to deal with the issue of the 
effectiveness of wall functions (in general terms) for separated flows. 

Ramaprian, Patel, and Choi (1978, 1980) measured surface pressures, friction 
coefficients, and velocity profiles for incompressible flow over a body of revolution 
at zero and ten degrees angle of attack. The body consisted of half of a prolate 
spheroid with a hemispherical nose cap. The zero degree angle of attack case was 
computed here using the k - c model with the present wall function formulation as 
well as the Baldwin- Lomax (1978) algebraic turbulence model and the Chien (1982) 
low-Reynolds-number k — e model. Each model was run using three different grids 
to demonstrate the advantage of wall functions for this type of flowfield. 

The experimental Reynolds number based on freestream velocity was 2 x 10 6 , 
and the freestream Mach number was approximately 0.06. The computations were 
run at the test Reynolds number, but the freestream Mach number was raised to 
0.10 to avoid the incompressible limit of the code while still maintaining essentially 
incompressible flow. In the experiment, the boundary layer was tripped near the nose 
(x = 0.04), so that the computations were run completely turbulent. 

Since the problem was symmetric about the centerline of the body, only half of 
the domain was solved. All grids were generated using the hyperbolic scheme of Chan 
and Steger (1991). Initially, an elliptic scheme was tried, but it proved impossible to 
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maintain the desired wall spacing. The hyperbolic scheme permitted excellent control 
over wall spacing, ran extremely fast, and was easy to use. 

All distances in this section are normalized by the length of the body. 


7.3.2 Fine grid 


For the first grid, a wall spacing was desired that would provide adequate resolu- 
tion for the Baldwin-Lomax and Chien models without requiring an excessive number 
of grid points. A wall spacing of A ~ 4 was chosen to provide at least one point in 
the viscous sublayer, giving a spacing at the wall of 5 x 10~ 5 in units of body length. 
The grid is stretched geometrically normal to the body using a stretching ratio of 
1-11. This relatively conservative value was chosen to minimize numerical errors due 


to grid stretching. 

Although the location of the outer edge of the domain cannot be precisely speci- 
fied when using hyperbolic grid generation, the desired value can be approximated. A 
value of twenty body lengths was chosen to minimize boundary proximity effects. A 
relationship between wall spacing, domain size, number of grid points, and geometric 
stretching ratio is given by Lewis and Pletcher (1986), 

1 )2 mai 


A z w = 


A'"- 1 - 1 


(7.11) 


where A z w is the wall spacing, z max is the domain size, n is the number of grid points, 
and K is the stretching ratio. Solving for n, 

_ Mi +&(*-»)] 

" = 1 + mTo 

For the present case, this equation yields n = 104. In the streamwise direction, 
101 grid points were deemed sufficient for reasonable resolution. The body in the 


(7.12) 
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experiment was suspended by wires and therefore had no sting, so an o-grid was 
the natural choice for the grid configuration. The outer boundary shape, which was 
determined by the hyperbolic grid generator, was close to a semicircle. The grid is 
shown in Figure 7.12, and the coordinate system is shown in Figure 7.13. At the 
leading edge x = 0 and at the trailing edge x — 1. z = 0 at the centerline of the 
body. 

At the outer edge of the domain, it is not clear (a priori) which points require 
an “inflow” boundary condition and which require an “outflow” boundary condition. 
This must be ascertained from the solution at each time step. At each boundary 
point, the inner product of the velocity vector and an outward normal vector to 
the boundary was computed. A value that is less than or equal to zero indicates an 
inflow point, and a value that is greater than zero indicates an outflow point. At inflow 
points, density and momentum components were fixed and pressure was extrapolated. 
At outflow points, pressure was fixed, and density and momentum components were 
extrapolated. At the axis, density and x-momentum were extrapolated. The z- 
momentum was set to zero due to symmetry. Flow along the axis streamline was 
irrotational, so Bernoulli’s equation was used to compute pressure. At the body 
surface, density and pressure were extrapolated, and momentum was set to zero. 

The procedure for testing for inflow or outflow was also used for the k — c 
equations. At inflow points both variables were specified, and at outflow points 
they were both extrapolated. At the axis, they were both also extrapolated. At the 
body surface the boundary conditions are part of the wall function formulation. For 
the Chien model, k and e are set to zero at the wall. Since the tunnel turbulence 
level was not measured in the experiment, a value of 1/2% ( k/U = 0.005) was 
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(a) Complete grid 



(b) Closeup of body 

Figure 7.12: Body of revolution, 101 x 104 grid 
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(c) Closeup of leading edge 



(d) Closeup of trailing edge 


Figure 7.12 (Continued) 



95 



Figure 7.13: Body of revolution coordinate system 


assumed. A freestream turbulent viscosity of 0.1 (nondimensional) was assumed, and 
freestream e was calculated from equation (3.6). 

Computed boundary-layer profiles are often affected adversely by even modest 
amounts of smoothing (Kaynak and Flores 1987). Since the physical shear stresses 
become large near walls, smoothing may be decreased in this region without causing 
the computation to become unstable. Various techniques of rolling off smoothing 
were tried, but the best solution was to simply turn the smoothing off at the seven 
grid points adjacent to the wall. 

The computed pressure coefficient distribution is compared with the experimen- 
tal data in Figure 7.14. All three turbulence models give good results, with a small 
discrepancy near x = 0.1. The discrepancy could be alleviated by increasing the grid 
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Figure 7.14: Pressure coefficient, body of revolution, 101 x 104 grid 


resolution in that region, but since it is small, the present grid was deemed adequate. 

Friction coefficients are shown in Figure 7.15. For the wall function case, the fric- 
tion coefficient is calculated directly from equation (4.1). For the other two models, 
equation (7.5) is used, with t w « fxdV/dn, where V is the velocity magnitude and n is 
the normal distance from the wall. The measurements from all three models compare 
well with the test data with the exception of the first measurement station, where the 
Baldwin- Lomax model gave values slightly below the others. The first measurement 
station is located in a region of adverse pressure gradient immediately after the flow 
accelerated over the nose in a strongly favorable pressure gradient. This is a difficult 
situation to simulate accurately, and, as will be seen in subsequent plots, none of the 
models reproduce the character of the velocity distribution very well at this point. 
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Figure 7.15: Friction coefficient, body of revolution, 101 x 104 grid 

The velocity distributions are shown in Figure 7.16. All models compare rea- 
sonably well with the test data. More importantly, the wall function solution is very 
close to the the solution from the other two models. It should be kept in mind that 
the grid point adjacent to the wall lies in or near the viscous sublayer, and not in the 
log region. In this respect, it may be thought of as a “worst case” test for the wall 
functions. 

7.3.3 Medium grid 

Another grid was generated with y+ « 16, with a wall spacing of 20 x 10 -5 . 
This was outside the optimum range for all three models, since the Baldwin-Lomax 
and Chien models should have points in the viscous sublayer, and wall functions are 
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Baldwin-Lomax 
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(a) x = 0.176 

(b) x = 0.241 
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U/Uinf 


(c) x = 0.333 (d) x = 0.426 

Figure 7.16: Velocity profiles, body of revolution, 101 x 104 grid 
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(a) x = 0.537 (b) x = 0.648 

Figure 7.16 (Continued) 

expected to work best when the point adjacent to the wall is in the log region. 

The ( grid lines are identical to those of the previous grid. The same domain size 
and stretching ratio were also used. Using equation (7.12), 91 grid lines are required 
in the normal direction. Closeups of the leading and trailing edge regions of this grid 
(101 x 91) are shown in Figure 7.17. 

Pressure coefficients are shown in Figure 7.18. They are nearly identical to those 
for the fine grid. 

Not all of the friction coefficients, shown in Figure 7.19, fare as well. The wall 
function solution still matches the test data well, but the other models cannot cope 
with this coarser wall spacing. 

The velocity profiles, shown in Figure 7.20, exhibit a similar trend, though not as 
pronounced. The Baldwin-Lomax and Chien solutions are reasonable, but are clearly 
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(a) Closeup of leading edge 



(b) Closeup of trailing edge 
Figure 7.17: Body of revolution, 101 x 91 grid 
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deficient as compared to the previous solution, while the wall function solution still 
looks good. 

7.3.4 Coarse grid 

A final grid was generated with Az+ « 48, or 60 x 10 -5 in physical coordinates. 
Closeups of the leading and trailing edge regions of this grid (101 x 80) are shown in 
Figure 7.21. The results for this grid are shown in Figures 7.22, 7.23, and 7.24. 

The 101 x 80 grid is clearly too coarse for the Chien and Baldwin-Lomax models, 
while the wall function solution still looks good. It should be noted that the wall 
function solutions for the different grids show some small differences. In the 101 x 80 
grid, the point adjacent to the wall is in the log region, while this is not true for 
the other two finer grids. The boundary conditions for k and e are more rigorous 
for points in the log region, as was discussed in the “Wall Functions” chapter. The 
different approaches are therefore expected to affect the results. 

7.4 Prolate Spheroid 


7.4.1 Introduction 

Computation of the flow over a prolate spheroid at angle of attack is a partic- 
ularly challenging problem for CFD. There exist regions of favorable and adverse 
pressure gradients, the flow on the leeward side may be massively separated, and 
laminar, transition, and turbulent regimes are frequently encountered. 

In the late 1970’s and early 1980’s, a series of experiments was carried out at 
DFVLR (Deutsche Forschungs- und Versuchsanstalt fur Luft- und Raumfahrt) to 
obtain detailed measurements of the surface flow on a prolate spheroid at angle of 
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(a) x = 0.537 (b) x = 0.648 


Figure 7.20 (Continued) 


attack (Kreplin, Meier, and Maier 1978; Meier and Kreplin 1980a; Meier and Kreplin 
1980b; Kreplin, Vollmers, and Meier 1982; Meier, Kreplin, and Vollmers 1983). One 
of the specific intents of the experiments was to provide data for turbulence modeling 
(Kreplin, Meier, and Maier 1978). Measurements were made over a range of Reynolds 
numbers and angles of attack. In some cases the boundary layer was tripped, while 
in others, natural transition was permitted to evolve. 

In addition to surface pressure measurements, hot films were applied to the 
surface of the body to directly measure the magnitude and direction of the wall shear 
stress. Since the hot film was quite thin, it may be assumed that it was located 
within the viscous sublayer. Knowledge of the velocity and viscosity at a specified 
point in the viscous sublayer may be used to deduce the shear stress. Each hot film 
probe contained two elements at different angles to the flow direction. The difference 
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(a) Closeup of leading edge 



(b) Closeup of trailing edge 
Figure 7.21: Body of revolution, 101 x 80 grid 
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Figure 7.22: Pressure coefficient, body of revolution, 101 x 80 grid 



Figure 7.23: Friction coefficient, body of revolution, 101 x 80 grid 
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U/Uinf 


(a) x = 0.176 


(b) x = 0.241 
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(c) x = 0.333 


U/Uinf 

(d) a: = 0.426 


Figure 7.24: Velocity profiles, body of revolution, 101 x 80 grid 
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(a) x = 0.537 



U/Uinf 

(b) x = 0.648 


Figure 7.24 (Continued) 


in response of each element was used to calculate the flow angle. 

A number of investigators have computed this flowfield using both boundary- 
layer methods and the Navier-Stokes equations, and compared the results with the 
DFVLR test data. The Navier-Stokes computations are summarized in Table 7.1. 

As can be seen in the table, a wide range of conditions were both tested and 
computed. The only computation shown which utilized a two-equation turbulence 
model was that of Kim and Patel. In their computation, a k — e model was employed 
with a one-equation model patched in near the wall. It was therefore of interest to 
carry out the present computation at one of the conditions computed by Kim and 
Patel in order to compare the performance of wall functions to a method utilizing a 
different wall treatment. 

The flow was laminar over most of the body in the Re Uao = 1.6 x 10 6 case, so 
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the case at a = 10° and Re Uoc = 7.2 x 10 6 was chosen for the present study. Another 
advantage of this case was that rather than tripping the boundary layer, it was 
permitted to go through natural transition to turbulence. This permitted trying a 
different method of forcing the transition point in the k — e model (as compared 
to Kim and Patel). A freestream Mach number of 0.20 was chosen to minimize 
compressibility effects without excessively hindering convergence. 

The use of wall functions is challenging for such a flow, due to the three- 
dimensionality of the boundary layer. In discussing their approach to solving this 
flowfield, Deng, Piquet, and Queutey (1990) stated that 

The wall function approach is avoided in this work so that the equations 
. . . are solved to the wall. For a significant increase in numerical troubles 
and of computing time (because the integration is carried out to y + = 

.1 — .3), the delicate problem of the threedimensional [stc] specification of 
the log-law is avoided. 

7.4.2 Grid 

It is desirable to have the grid point adjacent to the wall fall near the bottom 
of the log region. Since the values of Ay+ are not known a priori, the following 
method was used for an estimate. The definitions of y + and Cj (equations (4.3) 
and (7.5)) show that y + oc yjc] for a given freestream velocity and viscosity. Using 
the flat plate equation (7.4) for a rough estimate of friction coefficient, it may be 
shown that y + oc ln(-Re). In the coarse grid test case for the body of revolution, 
Re Uoo = 7.2 x 10 6 , the wall spacing was 60 x 10 -5 , and A y+ « 48. Using the 
above estimate of the variation of y + with Re Uoo , wall spacing of 60 x 10 -5 results 



Table 7.1: Previous computations of DFVLR prolate spheroid 
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in A y+ « 44 for the present Reynolds number. This value of Ay+ is in the desired 
range, so a wall spacing of 60 x 10 -5 is used for the present case. 

The test model was supported by a sting, and the sting is modeled in the present 
computation. A C-0 grid is therefore employed. The domain for the body of rev- 
olution test cases extended approximately 20 body lengths in each direction. This 
proved to be more than adequate, so a domain size of 10 body lengths was used for 
the present case. A stretching ratio of 1.15 was chosen normal to the wall, resulting 
in 57 grid points (from equation (7.12)). There are 121 grid points along the surface 
of the body and sting, and 53 in the circumferential direction. The circumferential 
lines are somewhat clustered toward the leeward side in order to improve resolution 
in the separated region. Due to symmetry, only half of the flowfield was computed. 
The grid, which was generated using the same hyperbolic method as described in the 
previous section, is shown in Figure 7.25. 

7.4.3 Boundary conditions 

The boundary conditions for the present case were the same as those used for the 
body of revolution in the previous section with the exception of the axis at the leading 
edge, and with the addition of the symmetry plane conditions. In the previous case, 
the flow was irrotational along the axis since the body was at zero angle of attack, 
and Bernoulli’s equation could be employed. In the present case this was not possible, 
so the following method was used. Density and the two components of momentum 
parallel to the symmetry plane were extrapolated to the axis. The component of 
momentum normal to the symmetry plane was set to zero due to symmetry. Using 
the extrapolated values of momentum, the pressure was then computed using the 
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(a) Entire domain, side view 



1 ( 11111111111 . 

(b) Closeup of body, side view 


Figure 7.25: Prolate spheroid, 121 x 53 x 57 grid 
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(c) Closeup of leading edge, side view 




(d) Closeup of trailing edge, side view 
Figure 7.25 (Continued) 
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(e) Entire domain, front view, near center of body 



(f) Closeup of body, front view, near center of body 
Figure 7.25 (Continued) 
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x-momentum equation. At the symmetry planes, values of all dependent variables 
are reflected, e.g. pk=\ = Pk=z- 

It was found in previous runs that, below a certain level, the freestream turbulent 
kinetic energy had little effect on the results. Larger values, however, had a marked 
effect. For the small values, any diffusion of turbulent kinetic energy in to the bound- 
ary is apparently overwhelmed by the production rate in the boundary layer. Since 
the “small” values gave the best agreement with the test data, A:/(uoo) 2 = 1 x 10 -8 
is used here. A small value is employed, rather than zero, because the source terms 
in the pe transport equation contain k in several denominators. A nondimensional 
freestream turbulent viscosity of 0.01 was assumed, and the freestream e was calcu- 
lated from equation (3.6). The small value of fi t was chosen so as not to affect the 
laminar flow significantly, as will be clarified below. As in the previous test cases, 
lower limiters on both k and e are set to freestream values. 

When using algebraic turbulence models, tripping a laminar boundary may be 
effected simply by turning the model on at a specified location in the flowfield. This 
is possible because the turbulent viscosity is a function of the local flowfield, and his- 
tory effects are not considered. With two-equation models, the turbulence quantities 
are transported throughout the domain. Kim and Patel (1991) solved the turbu- 
lence transport equations over the entire domain, and then, at the end of each step, 
overwrote the computed turbulent viscosities with zeros at those points designated 
as being laminar. A problem with this approach is that k and e are increasing in 
the boundary layer, mostly due to the production terms, in the region of laminar 
flow. Their profiles are quite developed by the time the points of forced transition 
are reached. In the present study, the time step was set to zero at all laminar grid 
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points. The freestream values of k and e were therefore maintained at these points, 
and the turbulence was permitted to evolve downstream of the trip line. Since a 
small freestream value of fj, t is used, the laminar flow was not significantly affected. 

The friction coefficient map taken from Kreplin, Vollmers, and Meier (1982) is 
shown in Figure 7.26. The vector lengths are proportional to the measured wall 
shear stresses, and they are oriented at the circumferential angles at which the shear 
stresses act. It is not possible to deduce a precise transition line from this map, but 
an approximate line may be determined from the regions where the friction coefficient 
rapidly increases. Figure 7.27 shows an unwrapped surface grid with points marked 
at the computational trip line. 

7.4.4 Additional considerations 

Smoothing coefficients of e 4 = 0.10 and e 2 = 0.028 were used (see equations (6.79) 
and (6.80)), based on several test runs to see how low they could be pushed without 
compromising stability. 

Friction coefficient angle data were measured in the experiment, based on the 
difference in the response of the two elements of the hot-film probes. The friction 
coefficient angle, which is assumed to be identical to the flow angle in the compu- 
tation, is changing as the wall is approached. The wall spacing is fairly large since 
wall functions are being used, and it is therefore difficult to deduce the limit of the 
angle at the wall. The flow angle was computed at the grid points near the wall, and 
second order Lagrangian extrapolation was used estimate the angle at the wall. 

A similar difficulty was encountered in computing the friction coefficient in the 
laminar region of the flowfield. Third order Lagrangian interpolation was used to 
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Figure 7.26: Friction coefficient map from Kreplin, Vollmers, and Meier (1982) 
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Figure 7.27: Unwrapped surface grid with trip points, prolate spheroid 
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estimate the velocity parallel to the wall, at a normal distance of 1 x 10 -5 . The 
shear stress was then computed from t w = nV p /n, where V p is the parallel velocity 
component, and n, the normal distance to the wall, is equal to 1 x 10 -5 . 

In Kreplin, Vollmers, and Meier (1982), both friction coefficients and friction 
coefficient angles were presented in the form of the friction coefficient map shown in 
Figure 7.26. It is difficult to deduce accurate friction coefficients and angles from data 
in this form. A tape containing the data had been sent to NASA Ames by Tuncer 
Cebeci of McDonnell-Douglas, and these are the data that were used for comparison. 
The tape also included surface pressures, which were not presented in the paper. 

7.4.5 Results 

Pressure coefficients at five axial locations are compared with the test data and 
with the computation of Kim and Patel (1991) in Figure 7.28. At all sections with 
the exception of the first, several test points fail to follow the trend of the remaining 
points. An example of this may be seen in Figure 7.28(d) at 9 « 115. None of 
these points coincide with the separation line, and no anomaly is apparent in the 
other measurements in the same regions, so the discrepancy is most likely due to 
measurement error. 

In the results from the present computation, a small kink may be observed in 
Figure 7.28(b) at 0 « 55 and in Figure 7.28(c) at 0 fa 25. These are the points at 
which the turbulence model was turned on in the computation to simulate transition. 
The pressure is affected by the sudden addition of the normal turbulent stresses (the 
— | Sijpk term in equation (3.1)). 

Both computations capture the character of the pressure distributions, but the 
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magnitudes of the pressures are somewhat off in both computations. Kim and Patel 
speculated that the differences were due to blockage effects in the wind tunnel or a 
bias in the pressure measurements. The present results imply that the differences 
may be due to errors in the computations, since the computed values straddle the 
measured values. 

Friction coefficients are shown in Figure 7.29. At all sections with the exception 
of the first, several test points fail to follow the trend of the remaining points. An 
example of this may be seen in Figure 7.29(d) at 9 « 115. None of these points coin- 
cide with the separation line, and no anomaly is apparent in the other measurements 
at these points. 

The match of the present computation is not very good at the first section. Here, 
the flow is fully laminar, and the difficulties in deducing laminar friction coefficients 
are encountered, as described above. This section is fairly close to the leading edge, 
and the boundary layer is relatively thin. The boundary layer is not well resolved 
here, encompassing approximately 5 — 7 grid points. 

The next section, at x = 0.395, has laminar flow near the windward symmetry 
plane and turbulent flow over the remainder of the circumference. The transition 
to turbulence after the k — e model is tripped is quite abrupt, even though the 
freestream k and e are maintained up to the trip point. The present computation 
shows the beginning of separation near 6 = 150°, but it is not yet evident in the test 
data. 

At x = 0.565, the present computation compares fairly well with the test data. 
The minimum C/, which is near the separation line, is close to that of the test data, 
and the increase in friction coefficient near the leeward symmetry plane is captured. 
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In the last two sections the dips in Cj are captured, but they are too close to the 
leeward symmetry plane. Neither computation does a very good job of predicting 
the magnitudes of the friction coefficients in the separated region. 

In three-dimensional separated flows, the definition of the separation line is not 
as straightforward as in two-dimensional flows, because the flow does not necessarily 
reverse in the streamwise direction. In the present flowfield, the separation line may 
be defined as the line where the circumferential velocity component changes direction. 
Examining Figure 7.26, it is interesting to note that the separation line does not 
correspond to the line of minimum friction coefficient. At each cross-section the 
separation point is most clearly shown in plots of friction coefficient angle, where the 
line crosses 7 = 0. The friction coefficient angle plots are shown in Figure 7.30. At the 




128 


first section, in which the flow is fully laminar (as may be seen in Figure 7.26), there is 
a separated region on the leeward side. At the second section, the flow on the leeward 
side has become turbulent, and it has reattached. As the flow moves downstream, it 
undergoes a second separation on the leeward side, and remains separated through 
the remaining measurement stations. The wall function computation captures the 
trend of the angles, but the magnitudes axe not accurate, particularly in the separated 
regions. Near the aft end of the spheroid, the computed separation line is too close 
to the leeward symmetry plane. In order to help clarify the location of the computed 
separation line, top and side views of computed and experimental oil flow patterns 
are shown in Figures 7.31 and 7.32. It is also evident in these plots that the computed 
separation line is too close to the leeward symmetry plane. 



















132 



(a) Computation 



(b) Experiment 


Figure 7.31: Surface oil flow pattern, top view 
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(a) Computation 



(b) Experiment 


Figure 7.32: Surface oil flow pattern, side view 
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8. CONCLUSIONS AND RECOMMENDATIONS 


A method has been developed for the application of wall functions to generalized, 
curvilinear, nonorthogonal grids. It has been applied to a series of test cases of vary- 
ing complexity. Flat plate test cases were first run to check for coding errors and to 
verify that the general formulation reduces correctly for this simple geometry. Com- 
puted friction coefficients were in good agreement with values from a semi-empirical 
equation. Since wall functions are based on the law of the wall, applicability of wall 
functions to these test cases was not an issue. 

The test cases for the body of revolution at zero angle of attack were more chal- 
lenging. This flowfield had favorable and adverse pressure gradients and curvature of 
the body surface. Pressure gradients were neglected in deriving the law of the wall, 
and if they are large, they may cause inaccuracies in the shear stress calculation. The 
standard k — e model is also known to be less accurate for flows with adverse pres- 
sure gradients and/or strong streamline curvature (e.g., De Henneau, Raithby, and 
Thompson 1990). Even with these limitations, it was shown that the wall function 
formulation gave results comparable to the Baldwin-Lomax algebraic model and the 
Chien low-Reynolds-number k — e model, all of which showed good agreement with 
test data when using a relatively fine grid. One of the primary advantages of wall 
functions is that they work well with coarse grid spacing, and this was also demon- 
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strated. As the grid spacing was increased, the wall functions continued to give good 
results while the other two models began to break down. The wall function solution 
did show a noticable (though not drastic) change as the grid spacing at the wall 
was increased. This was attributed to the different boundary conditions used for k 
and e in the viscous sublayer and the log region. Diminished grid resolution of the 
boundary layer may also have been a factor. 

Solution of flow over a prolate spheroid at angle of attack is an ambitious goal 
for any turbulence model, and particularly for wall functions. The use of the law 
of the wall is questionable for separated flows, since the friction velocity may not 
be an appropriate velocity scale. Also, it may be necessary to resolve the flow close 
to the wall to accurately capture the location of the separation line if the boundary 
layer is highly skewed. Even so, wall functions have been shown in many instances 
to work well for fairly complex flowfields, so their use for this case is not out of the 
question. The goal of the present study was to demonstrate a method for applying 
wall functions to general geometries. The effectiveness of wall functions for solving 
complex flowfields is a related but separate issue. 

The friction coefficients from the computation compared reasonably well with the 
test data in the regions of attached turbulent flow, and were comparable in accuracy 
to the two-equation model with a one-equation model patched in at the wall used 
by Kim and Patel (1991). Since the friction coefficient is computed directly from 
the wall function equations, it is the most important parameter to test the present 
formulation. The friction coefficient angle and separation line location did not match 
the test data as well as the friction coefficients, although the agreement was not 
unreasonable. 
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The problem of resolving the boundary layer in the laminar regions of the flow- 
field has severed potential solutions. One is simply to use a finer grid. One of the 
main reasons for using wall functions is to permit the use of coarser grids, thereby 
saving grid points, so this solution is not ideal. A grid could be chosen which is a 
compromise between the two requirements, however. Another solution would be to 
use varying wall spacing. The wall spacing could be finer near the nose of the body 
and gradually become coarser downstream. Alternatively, a separate grid with fine 
wall spacing could be patched into the nose region. The flowfield could then be solved 
using a grid-embedding technique such as chimera (Benek, Buning, and Steger 1985). 

There are several improvements which should be investigated to improve the wall 
function solution. It is possible to derive a modified law of the wall which accounts 
for streamwise pressure gradients (Ferrari 1959). The application of this equation to 
numerical schemes is not completely straightforward, since the equation takes on an 
indeterminate form as the pressure gradient approaches zero. There is a well defined 
limit at this point, however, so the problem should be surmountable. Although the 
pressure gradient is usually neglected in standard wall function formulations, its effect 
can be large (Mellor 1966; Tennekes and Lumley 1972), so the use of the modified law 
of the wall should improve the solution for flows with significant pressure gradients. 
It would be worthwhile to investigate this modification by solving a flowfield that is 
simpler than the prolate spheroid at angle of attack. A separating boundary layer in 
a two-dimensional diverging channel is a good choice. 

A related possible improvement to the present formulation involves the applica- 
tion of the computed shear stress to the right hand side of the Navier-Stokes equa- 
tions. The shear stress at a point between the wall and the grid point adjacent to 
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the wall is presently assumed to equal the wall shear stress. If the boundary layer 
assumption is invoked, the streamwise momentum equation evaluated at the wall 
yields the gradient of the shear stress as a function of streamwise pressure gradient. 
A better estimate of the shear stress near the wall is therefore possible. Although 
the boundary layer assumption is not correct in some regions of the flowfield, this 
modification should still be an improvement over neglecting the pressure gradient 
altogether. 

As mentioned above, the k — e model requires attention for flows with strong 
adverse pressure gradients and streamline curvature. Various modifications have been 
proposed (Launder, Priddin, and Sharma 1977; Hanjalic and Launder 1980; Rodi and 
Scheuerer 1983; Pourahmadi and Humphrey 1983). Evaluation of these modifications 
for relatively simple flowfields would be of value. The most promising formulation 
could then be applied to the prolate spheroid problem. The use of wall functions 
on curved surfaces should not pose any problems, since the law of the wall has been 
observed to hold close to both convex and concave surfaces (Moser and Moin 1987). 

Since neither the present computation nor the computation of Kim and Patel 
(1991) gave accurate friction coefficients in the separated region of the flowfield (with 
the exception of the middle measurement station in the present computation), the 
problem may be in the k — e model itself, rather than the wall treatment. Diagnosing 
this deficiency with the data available from the prolate spheroid measurements is 
difficult. It would be useful to have a detailed set of turbulence measurements (e.g. 
turbulent kinetic energy, Reynolds stress) available for a simpler separated flowfield 
such as flow over an airfoil or flow in a diverging duct. Computing the flowfield using 
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a variety of turbulence models would give insight into the specific deficiencies of each 
model. 
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10. APPENDIX A: FLUX JACOBIANS 


10.1 Navier-Stokes 


Before presenting these matrices, a word about their derivation is in order. Each 
flux vector may be considered to be a function of both the dependent variable vector 
Q and its spatial derivatives, for example dQjdx. We may therefore write for a 
one-dimensional equation 


A 


dE_ 

dQ 

dEjQM dEjQMdQ, 

dQ dQ x dQ 


( 10 . 1 ) 


It is common practice to neglect the second term on the right hand side (e.g. Pulliam 
1984), and this approximation will be made here. Using the notation of Pulliam, the 
Jacobian matrices are as follows. 


Msenm 


PRECEDING PAGE BLANK NOT FILMED 
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A = 


(t 6 

-u9 + £ x (j) 2 £t + 6-( 7 - 2 ) 

-u0 + £j,0 2 ^ x v-(7-1)4 w 

— + £ z <p 2 £ x w - (J - l) £ z u 
9[4> 2 ~a l \ ^xOl - (7 - 1)«0 


£*«“ ( 7 - l)£x” 
+ ^ - (7 - 2 ) £j,v 
ZyW-ii- l)£ z v 

Cy«l - (7 - l)v^ 


6 

0 


6« - ( 7 - i)^w 

(7-1)6. 


6*«-(7-l)6*«> 

(7- l)6y 

(10.2) 

£7 

4- 

«r> 

1 

1 

to 

s* 

8 

(7 - 1)6, 


- (7 - l)u>0 

7* + 6t 


where 



<N 

II 

<3 


(10.3) 

9 = + £yV + £ Z W 


(10.4) 

4> 2 = ^ (7 - 1 ) (w 2 + v 2 + w 2 ) 


(10.5) 

The flux Jacobians B and C are identical to the above matrix, with the 

exception 


that all occurrences of £ are replaced with rj or £ respectively. 
The viscous Jacobian, neglecting cross-derivative terms, is 


0 

0 

0 

0 

0 


m 2 1 

a, Sell 

Ul at 

o.fr 


0 


m 31 

^ dp-' 

Q 2ir 

a < S 6 T 

-Sr 

0 

( 10 . 6 ) 

m 41 

a ,^Ci 

“ 3 d( 

«»¥ 


0 


”*51 

m 52 

m 53 

77154 

”*55 
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where 



<*! — (/* + Vt) + fy + Cz) 


a 2 — « (A* "h Ah) 

03 = ^(A* + Ah)C*fz 
£*4 = (m + Ah) 

05 = |(m + /*»)?»6 


( 10 . 7 ) 

( 10 . 8 ) 

( 10 . 9 ) 


( 10 . 10 ) 

( 10 . 11 ) 

( 10 . 12 ) 

( 10 . 13 ) 

( 10 . 14 ) 

( 10 . 15 ) 

( 10 . 16 ) 

( 10 . 17 ) 

( 10 . 18 ) 

( 10 . 19 ) 

( 10 . 20 ) 
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«6 = in + Ht) (g +g + (10.21) 

As with the inviscid Jacobians, the viscous Jacobians for the other two coordinate 
directions may be obtained by substituting ij or ( for £ in the above equations. 


10.2 k-e 


For the standard k — e model, the source term Jacobian is 


where 


and 


2 e„A\ - B -c^) 2 - 1 

0^.4 + c 2 (f) 2 —c\B — 2c 2 ^ 

/ duj duj 2 duk\ dui 

dxi 3 u dxk } dxj 


B=- 


2 duk 

3 dxk 


( 10 . 22 ) 


(10.23) 


(10.24) 
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11. APPENDIX B: k - c SOLVER 


11.1 Banded Solver 


The k and e transport equations are coupled to one another through the diffusion 
terms and the source terms. In linearizing the diffusion terms, the turbulent viscosity 
was taken from the previous time step, so the only remaining coupling is through the 
source terms. The linearized source terms have been included in the £ direction 
factor, so the k and e equations in the f and r? direction factors are uncoupled from 
one another. A scalar banded solver may therefore be used for these factors. 

The form of the system of equations is 

Ax = B (11.1) 


where boldface print is used for matrices and vectors. The structure of the A matrix 
is shown in Figure 11.1, the vector of unknowns (x) in Figure 11.2, and the right 
hand side (B) in Figure 11.3. The lower diagonal is eliminated as follows: 


bki ~ bki c b.-i) j. 


K = “ %- 


q*. 

1} b 

%-u 

&k t 


B' ki = B k , - B^- 




(11.2) 

(11.3) 


(11.4) 
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h 2 0 
0 b f2 

a k3 0 

0 a e 3 


c *: 2 0 

0 c £3 

h 3 0 

0 b (3 

a k 4 0 

0 a f4 


c fc , 0 
0 c (3 


b k < 0 

0 b u 


c k4 0 
0 c (4 


Figure 11.1: Structure of banded matrix 


•^ki 

X (3 





Figure 11.2: Structure of vector of unknowns 
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This fully solves the system. 
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11.2 Block Solver 


The equations in the ( direction factor are coupled through the linearized source 
terms, so a block tridiagonal solver is required. Since the blocks are 2x2, they may 
be inverted algebraically. For a general nonsingular 2x2 matrix 


[Q] = 


q r 
s t 


the inverse is given by (e.g. Anton 1973) 


( 11 . 10 ) 


[ Q ]- 1 = 


qt — rs 


t —r 
-s q 


( 11 . 11 ) 


The procedure used here is a standard block tridiagonal solver with algebraic in- 
versions. The structure of the block matrix is shown in Figure 11.4. The solution 
procedure is similar to that for the banded scalar tridiagonal matrix. The lower 
diagonal is eliminated with the equations 


IK) = N - [a,] [b (l _i)] 1 [c^u] (11.12) 

and 

[B-] = [B.] - [a,] [b,,.,,]' 1 [B (i _„] (11.13) 

The backsweep begins in the last row, 

[x (A r — !)] = [b( Ar _ 1) ] [B^d] (11-14) 

For the remaining rows, 


[x,] = lb-]- 1 {[B-] - (cj [x„ +11 ]} 


(11.15) 
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[b 2 ] [c 2 ] 

N [b 3 ] [c 3 ] 

W [b 4 ] M 


Figure 11.4: Structure of block matrix 




